% 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 = 5000;
maxX = 100;
maxY = maxX;
f=0;
nNewParticlesPerFrame = 10;

% Initialize matrix containing all 2D gird points A(x,y)
% 1 <= x <= maxX
% 1 <= y <= maxY
% A(x,y)=0 ... cite is empty
% A(x,y)>0 ... cite is filled
A = ones(maxX, maxY);

% Introduce a solid wall of already aggreated particles at y=1
x0=maxX/2
y0=maxY/2
R=maxX/3;
for i=1:maxX;
    for j=1:maxY;
        dx=i-x0+1;
        dy=j-y0+1;
        r=hypot(dx,dy);
        if (r<R);
            A(i,j)=0;
        end
    end
end



% To save computer time, we want to inject the new particle not too far
% above growing aggreate. We inject at on a line 'yStart', which
% keeps being increased so that it is always 'yBuffer' lines above the
% highest frozen particle
yBuffer = 1;
yStart = 49+yBuffer;

for i = 1:nParticles;
    % Compute new starting point on line y=yStart
    x  = 50;
    y  = 50; %always start at upper limit
    
    while 1;
        xOld = x; %store current x and y values
        yOld = y;
        
        % Arbitrarily select direction to turn to
        r = rand();
        %%% Determine the direction of motion using the random number 'r'
        %%% Update 'x' and 'y' with the following probabilities
        %%% 25% go LEFT
        %%% 25% go RIGHT
        %%% 25% go UP
        %%% 25% go DOWN
        %%% --> need to enter about 10 lines of code
        
        if r<=0.25;
            y=y+1;
        end
        if 0.25<r && r<=0.5;
            x=x+1;
        end
        if 0.5<r && r<=0.75;
            y=y-1;
        end
        if 0.75<r && r<=1;
            x=x-1;
        end
        
        
        %%% Now apply periodic boundary conditions in X direction
        %%% x must stay with 1 <= x <= maxX
        %%% if x exits to the left, let it enter on the right
        %%% if x exits to the right, let it enter on the left
        %%% --> enter 6 lines of code (or less)
        
        
        
        if (A(x,y)>0);
            x=xOld;
            y=yOld;
            continue;
        end

        
        xR = x+1; % xRIGHT = x coordinate of my right neighbor
        xL = x-1; % xLEFT  = x coordinate of my left  neighbor
        yU = y+1; % yUP    = y coordinate of my upper neighor
        yD = y-1; % yDOWN  = y coordinate of my lower neighbo
        
      
        
        
        
        % Any of my 4 neighor filled already? If yes, get stuck here
        if (A(xR,y) > 0 | A(xL,y) > 0 | A(x,yU) > 0 | A(x,yD) > 0);
            p=rand();
            if p<0.1; %change sticking probability
                A(x,y) = 1; % mark site as take
                
                fprintf('i= %5d x=%4d y=%4d yStart=%4d\n',i,x,y,yStart);
                break;
            end
        end
        
    end
    
    if (i==nParticles | mod(i,nNewParticlesPerFrame)==0);
        imagesc(A');
        set(gca,'YDir','normal');
        colormap(1-gray);
        axis equal;
        axis tight;
        drawnow
        saveas(gcf, ['geode_point1_',sprintf('%03d',f),'.png'])
        f = f + 1;
    end
  
    if (yStart+1==maxY/2)
        fprintf('Growth reached Y limit after only %d particles\n',i);
        break;
    end
end

