% MATLAB code for drawing DLA by Chris Bishop
% Calls auxilary file plotcirc.m

% Draw DLA, but keep tract of the vertices of the convex
% hull at each stage. Draw DLA with two colors, depending
% on whether a new disk is on the convex hull boundary
% when it is added.
% This also draws the convex hull boundary of the DLA.

% This version plots multiple DLA clusters and collects
% some simple statistics on the rates of growth of the
% cluster and the vertices of its convex hull.

% n = number of disks to add to DLA
% N = number of different DLAs to draw

function dla_hull(n,N)
tic
saveK=zeros(N,n);
for iter=1:N
clear centers;
clear lograd;
tic;
centers=[0,2];
hullx=[0,2];hully=[0,0];
pointer=[0,0;1,2];
counter=2;
attempts=2;
% centers=zeros(1,10^7);
k=2;
while k < n
maxradius=max(abs(centers))+2;
z=(maxradius+1)* exp(2 * pi * i * random('unif',0,1));
step=min(abs(z-centers))-2;
while((step>.01)&&(step<5000))
z=z+ step * exp( 2 * pi * i * random('unif',0,1));
step=min(abs(z-centers))-2;
end % while
if (step > 4999) % use Poisson kernel to retart
end % if
if (step<.02)
k=k+1;
[m,index]=min(abs(z-centers));
pointer=[pointer;index,z];
attempts=attempts+1;
counter=counter+1;
centers(counter)=z;
newhullx=[hullx,real(z)];
newhully=[hully,imag(z)];
K=convhull(newhullx,newhully);
hullx=newhullx(K);
hully=newhully(K);
saveK(iter,counter)=length(K);
if max(K)==length(newhullx)
hull(counter)=1;
else
hull(counter)=0;
end % if
hull;
end % if
if mod(k,10000)==0
[k/10000,toc]
end % if
end % for k

% figure; hold on; axis equal;
% title('Red disks on convex hull when added');
% list=find(hull==0);
% L1=length(list);
% scatter(real(centers(list)),imag(centers(list)),'.','k');
% list=find(hull==1);
% L2=length(list);
% percentred=L2/(L1+L2)
% scatter(real(centers(list)),imag(centers(list)),'.','r');
% plot dla
% plotcirc(centers,hull);
% figure; hold on; grid on;
% title('Number of points on convex hull at time t');
% plot(saveK)
% L=length(saveK);
% plot([1:L],2.2*log([1:L]),'k','LineWidth',2);
new_trial=[iter, toc]
end % for iter
size(saveK)
figure
hold on; grid on;
title(['Number of convex hull verticies \\ averaged over ',num2str(m),' trials']);
list=sum(saveK)/N;
plot([100:n],list(100:n));
plot([100:n],2.2*log([100:n]),'k','LineWidth',2);
figure
hold on; grid on;
title(['Log plot of number of convex hull verticies \\ averaged over ',num2str(m),' trials']);
plot(log([100:n]),list(100:n));
save('saveK.Sept3.mat','saveK')
return