clear all
clc

% 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 = 30000;
maxX = 300;
maxY = maxX;

nNewParticlesPerFrame = 20;

% 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);
B = A;
% Introduce a solid wall of already aggregated 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;
particle_numb = 1;
frame = 0;
for i = 1:nParticles
	% Compute new starting point on line y=yStart
	x  = 1+floor(maxX*rand); 
	y = yStart; %always start at upper limit
%     y  = 1 + floor(maxY*yStart);
    
    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 && r <= .25
            x = xOld - 1;
        elseif r > .25 && r <= .5
            x = xOld + 1;
        elseif r > .5 && r <= .75
            y = yOld + 1;
        else
            y = yOld - 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 x > maxX
            x = 1;
        end
        if x < 1
            x = maxX;
        end
        
        if (y>yStart)   % did the particle diffuse through the ceiling?
            y = yStart; % just set back to max if it gets too high
        end
        
        if (A(x,y)>0) % ok, this site has been taken already, reject the move
            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
        
        %%% Now apply periodic boundary conditions to xL and xR
 
        if xR > maxX
            xR = 1;
        end
        if xL < 1
            xL = maxX;
        end
        
 		% Any of my 4 neighor filled already? If yes, count neighbors, get stuck here
        p = .03;
        
                 
        if (A(xR,y) > 0 || A(xL,y) > 0 || A(x,yU) > 0 || A(x,yD) > 0) 
                %         count number of neighbors filled n
            n = 0;
            count = 0;
            ran_prob = rand();
            if A(xR,y) > 0
                n = n + 1;
                count = count + B(xR,y);
            end
            if A(xL,y) > 0
                n = n + 1;
                count = count + B(xL,y);
            end
            if A(x,yU) > 0
                n = n + 1;
                count = count+ B(x,yU);
            end
            if A(x,yD) > 0
                n = n + 1;
                count = count + B(x,yD);
            end
            delay = particle_numb - (count/n);
            p = 1 - (2000 * delay)/nParticles;
            if ran_prob <= p
                A(x,y) = 1; % mark site as taken
                B(x,y) = particle_numb;
                particle_numb = particle_numb + 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
    end

%     if (i==nParticles || mod(i,nNewParticlesPerFrame)==0 || yStart+1==maxY) 
%         figure(1)
%         imagesc(A');
%         set(gca,'YDir','normal');
%         colormap(1-gray);
%         axis equal;
%         axis tight;
%         drawnow
%     end
    
    if (i==nParticles || mod(i,nNewParticlesPerFrame)==0 || yStart+1==maxY) 
        frame = frame + 1;
        figure(2)
        imagesc(B');
        set(gca,'YDir','normal');
%         colormap(winter);
        axis equal;
        axis tight;
        drawnow
        saveas(gcf, ['movie_1_frame_',sprintf('%03d',frame),'.png']) 
    end

    if (yStart+1==maxY) 
        fprintf('Growth reached Y limit after only %d particles\n',i);
        break;
    end
end

	