% DIFFUSION LIMITED AGGREGATION along a sticky wall using Brownian motion
% 
% Read "The Science of Fractal Images", Ed. Peitgen and Saupe, p. 37 (1988)

% Author : B. Militzer, University of California, Berkeley 
% Date   : Sep 21, 2009

nParticles = 10000;

nNewParticlesPerFrame = 100;

n = 0;

t = linspace(-pi,pi, 350);
X = t .* sin( pi * sin(t)./t);
Y = -abs(t) .* cos( pi * sin(t)./t);

A = zeros(length(t));
maxX = size(A,1);
maxY = size(A,2);
centerx = maxX/2;
centery = maxY/2 + 50;

disty = linspace(-3.5,1,length(Y));
distx = linspace(-2.5,2.5,length(X));
spacey = disty(2)-disty(1);
spacex = distx(2)-distx(1);

for i = 1:length(t)
    xfill = X(i);
    yfill = Y(i);
    xind = find(distx < xfill + 4*spacex & distx > xfill - 4*spacex);
    yind = find(disty < yfill + 4*spacey & disty > yfill - 4*spacey);
    A(xind,yind) = 100; 
end

yBuffer = 5;
yStart = 130+yBuffer;
count = 0;
for i = 1:nParticles
    x = centerx;
    y = centery;
    while 1
        if A(centerx,centery) > 0
            fprintf('Growth reached center limit after only %d particles\n',i);
            break;
        end
        xOld = x; 
        yOld = y;
        
        r = rand();
     
        if 0<=r && r<0.3; 
            x = x + 1;
        elseif 0.3<=r && r<.6; 
            x = x - 1;
        elseif 0.6<=r && r<0.8; 
            y = y + 1;
        elseif 0.8<=r && r<=1; 
            y = y - 1;
        end
   
    if x>maxX; x=1; end
    if x<1; x=maxX; end
  
        
if (y>yStart)   
    y = yStart; 
end
if y < 1; y=1; end 
if (A(x,y)>0) 
    x=xOld;
    y=yOld;
    continue;
end

xR = x+1; 
xL = x-1; 
yU = y+1; 
yD = y-1; 
if xR>maxX; xR=1; end
if xL<1; xL=maxX; end
if yD < 1; yD = 1; end 

pstick = 1; q = 0.5;
if (A(xR,y) > 0 || A(xL,y) > 0 || A(x,yU) > 0 || A(x,yD) > 0) && q < pstick;
    count = count+1;
    A(x,y) = mod(count,500); 
    if (y+yBuffer>yStart && y+yBuffer<maxY)
        yStart = y+yBuffer;
    end
    fprintf('i= %5d x=%4d y=%4d yStart=%4d\n',i,x,y,yStart);
    break;
end

    end
    
    if (i==nParticles || mod(i,nNewParticlesPerFrame)==0) %|| yStart+1==maxY)
        imagesc(A');
        set(gca,'YDir','normal');
        colormap('hot')
        caxis([0 500])
        axis equal;
        axis tight;
        drawnow
        n = n+1;
        filename = ['2ndheart_frame_number_',sprintf('%03d',n),'.png']
        print(filename,'-r600','-dpng')
    end
    
    % if A(centerx,centery) > 0
    %    fprintf('Growth reached center limit after only %d particles\n',i);
    %   break;
    % end
end
title(['Particle only sticks with p = ' num2str(pstick)])

