% 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 = 1000000;
maxX = 100;
maxY = 400;
maxZ = 100;

nNewParticlesPerFrame = 100; %%makes code runs faster

% 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,maxZ);   %%sets matrix to all zeros
B = zeros(maxX, maxZ);
% Introduce a solid wall of already aggreated particles at y=1 
A(:,1,:) = 1;  %%set first row to 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 %%y start is a point where the particles are
% introduced
yBuffer = 5;  %%5 pixles over the trees
yStart = 100; %%%where is starts

for i = 1:nParticles
	% Compute new starting point on line y=yStart
	x  = 1+floor(maxX*rand); %%rand numbers between 0 and 1, set to an integer, floor rounds the integer down
	y  = yStart; %always start at upper limit
	z = 1+floor(maxZ*rand);

    while 1  %%infinite loop, needs a break
        xOld = x; %store current x and y values
        yOld = y;
	    zOld = z;
        
        % Arbitrarily select direction to turn to
		r = rand();
        if r<0.025;
            x=x-1;
        elseif r<=0.05;
             x=x+1;
        elseif r<0.075;
             y=y+1;
	elseif r<0.1;
	     z=z+1;
	elseif r<0.125;
	     z=z-1;
	elseif r<0.15;
	     z=z+1;
	     x=x+1;
	elseif r<0.175;
	     z=z+1;
	     x=x-1;
	elseif r<0.2;
	     z=z-1;
	     x=x+1;
	elseif r<0.225;
	     z=z-1;
	     x=x-1;
	elseif r<0.24;
	     y=y-1;
	     z=z+1;
	     x=x+1;
	elseif r<0.32;
	     y=y-1;
	     z=z+1;
	     x=x-1;
	elseif r<0.40;
	     y=y-1;
	     z=z-1;
	     x=x+1;
	elseif r<0.48;
	     y=y-1;
	     z=z-1;
	     x=x-1;
	elseif r<0.56;
	     y=y-1;
	     z=z+1;
	elseif r<0.64;
	     y=y-1;
	     z=z-1;
	elseif r<0.72;
	     y=y-1;
	     x=x+1;
	elseif r<0.8;
	     y=y-1;
	     x=x-1;
        else
             y=y-1;
        end
        %%if statements
        %%% 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
        
        ..................................................
             
        %%% 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 it exceeds maxX/minX
        %%% 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 z>maxZ
	    z=1;
	end
	if z<1
	    z=maxZ;
	end
        
        

        ..................................................
        
        if (y>yStart)   % did the particle diffuse through the ceiling?
            y = yStart; % just set back to max if is gets too high
        end
        
        if (A(x,y,z)>0) % ok, this site has been taken already, reject the move
            x=xOld;
            y=yOld;
	    z=zOld;
            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
	zB = z+1;
	zT = z-1;
        
        %%% Now apply periodic boundary conditions to xL and xR
 
        if xR>maxX
            xR=1;
        end
        if xL<1
            xL=maxX;
        end
	if zB>maxZ
	    zB=1;
	end
	if zT<1
	    zT=maxZ;
	end
            
        
        ..................................................
        
 		% Any of my 4 neighor filled already? If yes, get stuck here
        if (A(xR,y,z) > 0 | A(xL,y,z) > 0 | A(x,yU,z) > 0 | A(x,yD,z) > 0 | A(x,y,zB) > 0 | A(x,y,zT) > 0) %#ok<*OR2>
            A(x,y,z) = 1; % mark site as taken

	if A(x,y,z)==0
            p=0.84;
            r=rand();
            if p>r
                A(x,y,z)=1;
            end
        end
            if (y+yBuffer>yStart && y+yBuffer<maxY) 
                yStart = y+yBuffer;
            end
            %fprintf('**i= %5d x=%4d y=%4d z=%4d yStart=%4d\n',i,x,y,z,yStart);
            break;
        end

    end

    if (i==nParticles | mod(i,nNewParticlesPerFrame)==0 | yStart+1==maxY)
        fprintf('i = %4d\n',i);
        for j = 1:maxX;
            for k = 1:yStart;
                for h = 1:maxZ;
                    B(j,k) = (B(j,k)*maxZ + A(j,k,h)) / maxZ;
                end
            end
        end
        imagesc(B');
        set(gca,'YDir','normal');
        colormap(bone);
        axis equal;
        axis tight;
        xlabel('Width','FontSize',10)
        ylabel('Height','FontSize',10)
        set(get(gca,'ylabel'),'Rotation',0.0)
        title('3D DLA','FontSize',14)
        drawnow
    end

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

	