-
Notifications
You must be signed in to change notification settings - Fork 0
/
DLARectangular2.m
196 lines (159 loc) · 5 KB
/
DLARectangular2.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
%% DLA SIMULATION USING SQUARE MATRIX
clear
clf
clc
% Set up variables
width = input('width: ');
tic
% numberOfParticles = input('Number of Particles: ');
%% Initial Setup:
% We create a matrix of zeros and then plant a seed at the centre of the
% matrix. Change entry of matrix to 1 if it is aggregated.
%
% This is the matrix set up, we start from inner square and consider
% particle escaped if it reaches the outer square
% -----------
% | ---- |
% | | | |
% | | | |
% | ---- |
% -----------
matrix = zeros(2*width,2*width);
middlex = width;
middley = width;
matrix(middley,middlex) = 0.1;
% We want to let a particle randomly enter from the edges of the matrix so
% first let's choose a side
particleNumber = 1; %integer
x = 0; %integer
y = 0; %integer
aggregate = 0; %bool
endScript = 0; %bool
escape = 0; %bool
%% Random Walk Script
while endScript == 0
%first decide where the particle starts from
randSide = rand;
if rand <0.25
%start on bottom
randPlace = rand*width + width/2;
randPlace = ceil(randPlace);
x = randPlace;
y = ceil(3*width/2);
elseif rand <0.5
%start on right
randPlace = rand*width + width/2;
randPlace = ceil(randPlace);
x = ceil(3*width/2);
y = randPlace;
elseif rand<0.75
%start on top
randPlace = rand*width + width/2;
randPlace = ceil(randPlace);
x = randPlace;
y = ceil(width/2);
else
%start on left
randPlace = rand*width + width/2;
randPlace = ceil(randPlace);
x = ceil(width/2);
y = randPlace;
end
% Now we have a starting point on the edge of the inner square of the matrix. Now we want
% to perform random walk in the matrix until we are next to an entry
% which is 1. Set aggregate to false.
aggregate = 0;
escape = 0;
while (aggregate == 0) && (escape == 0)
randWalk = rand;
if randWalk < 1/4
% go right
x = x+1;
elseif randWalk < 2/4
% go left
x = x-1;
elseif randWalk < 3/4
% go up
y = y-1;
else
% go down
y = y+1;
end
% escape test
if (x == 1) || (x == 2*width) || (y == 1) || (y == 2*width)
escape = 1;
end
% aggregate test, only if not escaped
if escape == 0
if (matrix(y,x+1) + matrix(y,x-1) + matrix(y+1,x) + matrix(y-1,x)) ~= 0
aggregate = 1;
end
end
%We keep repeating this until the particle is finally stuck
end
if escape
disp('particle escaped')
end
if aggregate
particleNumber = particleNumber + 1;
disp(['particle aggregated: ' num2str(particleNumber)])
if particleNumber < 500
matrix(y,x) = 0.1;
elseif particleNumber < 1000
matrix(y,x) = 0.2;
elseif particleNumber < 1500
matrix(y,x) = 0.3;
elseif particleNumber < 2000
matrix(y,x) = 0.4;
elseif particleNumber < 2500
matrix(y,x) = 0.5;
elseif particleNumber < 3000
matrix(y,x) = 0.6;
elseif particleNumber < 3500
matrix(y,x) = 0.7;
elseif particleNumber < 4000
matrix(y,x) = 0.8;
elseif particleNumber < 4500
matrix(y,x) = 0.9;
else
matrix(y,x) = 1.0;
end
% end of script check
if (x==ceil(width/2) || x==ceil(3*width/2) || y==ceil(3*width/2) || y==ceil(width/2))
endScript = 1;
% now want to calculate diameter of fractal
% we could change to make sure that we always add particles until we hit
% the edge, then diameter we could work out
xdif = abs(x - middlex);
ydif = abs(y - middley);
radius = sqrt(ydif^2 + xdif^2);
diameter = 2*radius;
fractdim = log(particleNumber)/log(radius);
end
end
end
%Tell why escaped, if true
% if (endScript == 1)
% disp ('Script ended early as aggregation was about to exceed matrix dimensions');
%
% end
% uncomment this for diagram to just appear at end
%% Plot graph
figure(1)
imagesc(matrix)
colormap(jet)
title(['DLA with ' num2str(particleNumber) ' particles'])
%text(width/3,3*R/2 + R/6,['Fractal Dimension: ' num2str(fractalDimension)]);
%text(R/3,3*R/2,['Radius: ' num2str(maximumDistance)]);
timeElapsed = toc;
%text(R/3,3*R/2 - R/6, ['Time Elapsed: ' num2str(timeElapsed)]);
%xlabel(num2str(2*width))
%ylabel(num2str(2*width))
axis equal
xlim([floor(width/2),ceil(3*width/2)])
ylim([floor(width/2),ceil(3*width/2)])
%% Display Outputs
disp(['Number of particles: ' num2str(particleNumber)]);
disp(['Fractal Dimension: ' num2str(fractdim)]);
disp(['Diameter: ' num2str(diameter)]);
disp(['Time Elapsed: ' num2str(timeElapsed)]);