function diffusion_limited_aggregation2D(nParticles, maxX, maxY)

% DIFFUSION LIMITED AGGREGATION along a wall using Brownian motion
% 
% Examples:
%    diffusion_limited_aggregation2D(1000,200,50) 
%    diffusion_limited_aggregation2D(3000,200,100) 
%    diffusion_limited_aggregation2D(10000,500,200) 

% Read "The Science of Fractal Images", Ed. Peitgen and Saupe, p. 37 (1988)

% Author : B. Militzer, University of California, Berkeley 
% Date   : Aug 25, 2008
	
if nargin < 3
	error('Type "diffusion_limited_aggregation2D" for calling syntax.');
end
	
% 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 = zeros(maxX, maxY);
% Introduce a solid wall of already aggreated particles at y=1 
A(:,1) = 1;

% 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 = 5;
yStart = 1+yBuffer;

for i = 1:nParticles
	% Compute new starting point on line y=yStart
	x  = 1+floor(maxX*rand); 
	y  = yStart; %always start at upper limit

    while 1
        % Arbitrarily select direction to turn to
		rval = rand;
        %%% Determine the direction of motion using the random nubmer 'rval'
        %%% 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
                
        .....................................................
            
        %%% 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 y exits to the right, let it enter on the left
        %%% --> enter 6 lines of code (but is also possible with one using 'mod' function) 

        .......................................................
        
        if (y>yStart)
            y = yStart; % just set back to max if is gets too high
        end
               
        xp = x+1;
        xm = x-1;
        xp = mod(xp-1,maxX)+1; % efficient way to apply p.b.c.
        xm = mod(xm-1,maxX)+1; % efficient way to apply p.b.c.
        
 		% Any frozen particle nearby? If so attach here 
        if (A(xp,y) > 0 | A(xm,y) > 0 | A(x,y+1) > 0 | A(x,y-1) > 0)
            A(x,y) = 1;
            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 (yStart+1==maxY) 
        fprintf('Growth reached Y limit after only %d particles\n',i);
        break;
    end
end

imagesc(maxX:1,1+maxY-(1:maxY),A');
colormap(1-gray);
axis equal;
axis tight;