% 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

clear
clf

%Make figure full screen size to maximize video resolution
scrsz = get(0,'ScreenSize');
figure('Position',[0 0 scrsz(3) scrsz(4)-50])

set(gca,'nextplot','replacechildren');

aviobj = avifile('sandy_movie_final','quality',100,'compression','Cinepak');
q=100;

maxX = 500; 
maxY = 500;

xstart=125;
xend=375;
ystart=xstart;
yend=xend;

nParticles = ((xend-xstart)+1)*((yend-ystart)+1); 

nNewParticlesPerFrame = 100;
 
A = zeros(maxX,maxY);
A(xstart:xend,ystart:yend) = 1;

while nParticles>0;
    for y=ystart:yend;
        if mod(y,2)==0;
            dx=0;
        else dx=1;
        end
        for x=xstart+dx:2:xend;
         ..................................................
if A(x,y)>0;

        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 neighbor
   
% % %     Could have also written the part below in a much easier way: if A(xR,y)>0, n=1, end; if A(xL,y)>0, n=n+1, end; etc.     
% % %     I am trying to say A(x,y)=0 with r probability.
    r=rand();
    if (A(xR,y)<1 && A(xL,y)<1 && A(x,yU)<1 && A(x,yD)<1)&& r<=1; %If i put the probabilities in the if statement, can I take out the A combinations.
        %n=4; 4 neighbors have dissolved;
        A(x,y)=0;
        nParticles=nParticles-1;
        continue
    end
    
    if (A(xR,y)<1 && A(xL,y)<1 && A(x,yU)<1) || (A(xR,y)<1 && A(xL,y)<1 && A(x,yD)<1) || (A(xR,y)<1 && A(x,yU)<1 && A(x,yD)<1) || (A(xL,y)<1 && A(x,yU)<1 && A(x,yD)<1)&& r<=0.75;
        %n=3;
        A(x,y)=0;
        nParticles=nParticles-1;
        continue
    end
    
    if (A(xR,y)<1 && A(xL,y)<1) || (A(x,yU)<1 && A(x,yD)<1) || (A(xR,y)<1 && A(x,yU)<1) || (A(xR,y)<1 && A(x,yD)<1) || (A(xL,y)<1 && A(x,yU)<1) || (A(xL,y)<1 && A(x,yD)<1) && r<=0.5;
        %n=2;
       A(x,y)=0;
       nParticles=nParticles-1;
        continue
    end
    
    if (A(xR,y)<1 ||A(xL,y)<1 || A(x,yU)<1 || A(x,yD)<1) &&  r<=0.25;
        %n=1;
        A(x,y)=0;
        nParticles=nParticles-1;
        continue
    end
end
        end
    end
  
    for y=ystart:yend;
        if mod(y,2)==1;
            dx=0;
        else dx=1;
        end
        for x=xstart+dx:2:xend;
         ..................................................
if A(x,y)>0;

        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 neighbor
   
% % %     Could have also written the part below in a much easier way: if A(xR,y)>0, n=1, end; if A(xL,y)>0, n=n+1, end; etc.     
    r=rand();
    if (A(xR,y)<1 && A(xL,y)<1 && A(x,yU)<1 && A(x,yD)<1)&& r<=1; %If i put the probabilities in the if statement, can I take out the A combinations.
        %n=4; 4 neighbors have dissolved;
        A(x,y)=0;
        nParticles=nParticles-1;
        continue
    end
    
    if (A(xR,y)<1 && A(xL,y)<1 && A(x,yU)<1) || (A(xR,y)<1 && A(xL,y)<1 && A(x,yD)<1) || (A(xR,y)<1 && A(x,yU)<1 && A(x,yD)<1) || (A(xL,y)<1 && A(x,yU)<1 && A(x,yD)<1)&& r<=0.75;
        %n=3;
        A(x,y)=0;
        nParticles=nParticles-1;
        continue
    end
    
    if (A(xR,y)<1 && A(xL,y)<1) || (A(x,yU)<1 && A(x,yD)<1) || (A(xR,y)<1 && A(x,yU)<1) || (A(xR,y)<1 && A(x,yD)<1) || (A(xL,y)<1 && A(x,yU)<1) || (A(xL,y)<1 && A(x,yD)<1) && r<=0.5;
        %n=2;
       A(x,y)=0;
       nParticles=nParticles-1;
        continue
    end
    
    if (A(xR,y)<1 ||A(xL,y)<1 || A(x,yU)<1 || A(x,yD)<1) &&  r<=0.25;
        %n=1;
        A(x,y)=0;
        nParticles=nParticles-1;
        continue
    end
end
        end
    end

    q=q+1;

        imagesc(A);
        set(gca,'YDir','normal');
        colormap(1-gray);
        axis equal;
        axis tight;
        drawnow
    
    if mod(q,1)==0;
    fprintf('Frame %4d of %4d\n',q,100);
    F(q) = getframe;
    aviobj = addframe(aviobj,F(q));
    end
     
display(nParticles)    
end
aviobj = close(aviobj);
     
    
