nParticles = 10000;
maxX = 200;
maxY = 300;

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);

%sugar!
sugarX=maxX/2;
sugarY=5*maxY/6;
for i=(2*maxX/5):(3*maxX/5)
    for j=(2*maxY/3):maxY
        if abs(i-sugarX)+abs(j-sugarY) < maxX/12;
            A(i,j)=-2;
        end
    end
end

% Make the absinthe glass and liquid area.  Sorta.
for i=1:maxX
    for j=1:(2*maxY/3)
        if j < 2*maxY/3 - 10
            A(i,j)=-1;
        end
        if j < (((i-(maxX/2))/7)^2 + 50)
            A(i,j)=1;
        end
    end
end

% 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

yStart = round(2*maxY/3) - 10;

for i = 1:nParticles
	% Compute new starting point on line y=yStart
    if mod(i,2)==0
        x  = round(2*maxX/5);
    else
        x  = round(3*maxX/5);
    end
	y  = yStart; %always start at upper limit

    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 0<=r && r<=0.25;
            x=xOld+1;
            %y=yOld
        end
        if 0.25<=r && r<=0.5;
            x=xOld-1;
            %y=yOld
        end
        if 0.5<=r && r<=0.73;
            %x=xOld
            y=yOld+1;
        end
        if 0.73<=r && r<=1;
            %x=xOld
            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 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 xR>maxX;
            xR=1;
        end
        if xL<1;
            xL=maxX;
        end
        
 		% Any of my 4 neighor filled already? If yes, get stuck here
        if (A(xR,y) > 0 | A(xL,y) > 0 | A(x,yU) > 0 | A(x,yD) > 0)
            A(x,y) = 1; % mark site as taken
            %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(1-gray);
        axis equal;
        axis tight;
        drawnow
        saveas(gcf, ['absinthe_',sprintf('%03d',i/nNewParticlesPerFrame),'.png'])
    end

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

	