clear;
clf;

% 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 = 20000;
maxX = 70;
maxY = 70;

nNewParticlesPerFrame = 100;

% 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:2:maxX,1) = 1;
A(2:2:maxX,1) = 2;

yBuffer = 5;
yStart = 1+yBuffer;

type=1


for i = 1:nParticles
	% Compute new starting point on line y=yStart
	x  = 1+floor(maxX*rand); 
	y  = yStart; %always start at upper limit
othertype=3-type;
    
    while 1
        xOld = x; %store current x and y values
        yOld = y;
        
		r = rand();

        if 0<r & r<0.25;
           x=x-1;
        end;
        if 0.25<r & r<0.50;
           x=x+1;
        end;
        if 0.50<r & r<0.75;
           y=y-1;
        end;
        if 0.75<r & r<1;
           y=y+1;
        end;
        
      
             
        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 is 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 xL<1;
           xL=maxX;
        end;
        if xR>maxX;
           xR=1;
        end;
        
 %no=n othertype
 %nt=n type
        no=0;
        nt=0;
       if A(xR,y)==othertype;
          nt=no+1;
       end
       if A(xL,y)==othertype;
          no=no+1;
       end
       if A(x,yD)==othertype;
          no=no+1;
       end;
       if A(x,yU)==othertype;
          no=no+1;
       end;
       if A(xR,y)==type;
          nt=nt+1;
       end
       if A(xL,y)==type;
          nt=nt+1;
       end
       if A(x,yD)==type;
          nt=nt+1;
       end;
       if A(x,yU)==type;
          nt=nt+1;
       end;
       
       % Probabilities
       p=0.03;
%       if no==0 && (nt==3 || nt==4);
       if (no==0) 
          p=0;
       end    
       if no==0 && nt==2;
          p=10*p;
       end
       if no==0 && nt==1;
          p=10*p;
       end
       if no==1 && nt==3;
          p=10*p;
       end
       if no==1 && nt==2;
          p=5*p;
       end
       if no==1 && nt==1;
          p=p;
       end
       if no==1 && nt==0;
          p=p;
       end
       if no==2 && nt==2;
          p=2*p;
       end
       if no==2 && nt==1;
          p=p;
       end
       if no==2 && nt==0;
          p=4*p;
       end
       if no==3 && nt==1;
          p=p;
       end
       if (no==3 || no==4)&& nt==0;
          p=p;
       end
%       fprintf('no=%d nt=%d p=%f\n',no,nt,p);
%       pause
       
        r2=rand();
        if r2<p;
            
            A(x,y) = type;
            type=othertype; % mark site as taken
           
            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);
        axis equal;
        axis tight;
        drawnow
    % end

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

end

	