clear all
clc
% 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 = 30000;
maxX = 300;
maxY = maxX;
nNewParticlesPerFrame = 20;
% 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);
B = A;
% Introduce a solid wall of already aggregated 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;
particle_numb = 1;
frame = 0;
for i = 1:nParticles
% Compute new starting point on line y=yStart
x = 1+floor(maxX*rand);
y = yStart; %always start at upper limit
% y = 1 + floor(maxY*yStart);
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 r > 0 && r <= .25
x = xOld - 1;
elseif r > .25 && r <= .5
x = xOld + 1;
elseif r > .5 && r <= .75
y = yOld + 1;
else
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 it 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, count neighbors, get stuck here
p = .03;
if (A(xR,y) > 0 || A(xL,y) > 0 || A(x,yU) > 0 || A(x,yD) > 0)
% count number of neighbors filled n
n = 0;
count = 0;
ran_prob = rand();
if A(xR,y) > 0
n = n + 1;
count = count + B(xR,y);
end
if A(xL,y) > 0
n = n + 1;
count = count + B(xL,y);
end
if A(x,yU) > 0
n = n + 1;
count = count+ B(x,yU);
end
if A(x,yD) > 0
n = n + 1;
count = count + B(x,yD);
end
delay = particle_numb - (count/n);
p = 1 - (2000 * delay)/nParticles;
if ran_prob <= p
A(x,y) = 1; % mark site as taken
B(x,y) = particle_numb;
particle_numb = particle_numb + 1;
if (y+yBuffer>yStart && y+yBuffer