nParticles = 20000;
maxX = 400;
maxY = 200;
p = 1;
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) = 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;
for i = 1:nParticles
% Compute new starting point on line y=yStart
x = 1+floor(maxX*rand);
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();
if (r < 0.4);
x = x+1;
elseif (r < 0.8);
x = x-1;
elseif (r < 0.9);
y = y+1;
elseif (r < 1);
y = y-1;
end
%%% 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 x exits to the right, let it enter on the left
%%% --> enter 6 lines of code (or less)
if x <= 0;
x = maxX;
end
if x >= maxX;
x = 1;
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 == 0;
xL = maxX;
end
..................................................
s = rand();
n= A(xR,y)+A(xL,y)+A(x,yU)+A(x,yD);
% Any of my 4 neighor filled already? If yes, get stuck here
if (n>0)
p=1;
if n>1
p=10*(n-1)*p;
end
if syStart && y+yBuffer