anArray = zeros(1,10); %Initialize storage array
storageCounterMax = 10;
for storageCounter = 1:storageCounterMax
%Variable Declaration
n = 1; %Initialize loop counter
N=0.25; %Range from random selection will be [-N,N]
loopMax=10000; %Number of random walk movements
ColorSet = varycolor(loopMax); %Set amount of colours
diameter = 6.94; %Set diameter of cylinder
radius = diameter*0.5; %Set radius of cylinder
totalHeight = 10.85; %Set total height of container
pauseTime = 0.0; %Set time between plotting iterations
boundaryCondition = true;
%Cylinder Figure
maxHeight = totalHeight - radius;
h=maxHeight;
x0=0;y0=0;z0=0;
[x,y,z]=cylinder(radius);
x=x+x0;
y=y+y0;
z=z*h+z0;
%Half Sphere Figure
A = [0 0 0 radius];
[X, Y, Z] = sphere;
XX = X * A(4) + A(1);
YY = Y * A(4) + A(2);
ZZ = Z * A(4) + A(3);
XX = XX(11:end,:);
YY = YY(11:end,:);
ZZ = -ZZ(11:end,:);
%Display Information
figure
xlabel('X');
ylabel('Y');
zlabel('Z');
surface(x,y,z, 'FaceAlpha', 0.1, 'FaceColor', [ 1 1 0], 'EdgeColor', [0.4, 0.4, 0.4]); hold on
surface(XX, YY, ZZ, 'FaceAlpha', 0.1, 'FaceColor', [ 1 1 0], 'EdgeColor', [0.4, 0.4, 0.4]); hold on
axis equal
view(3);
set(gca,'GridLineStyle','-')
Information=strcat('Brownian Motion in 3D Space');
title(Information ,'FontWeight','bold');
view(-109,40);
hold all
set(gca, 'ColorOrder', ColorSet);
%Initialize Beginning Location
heightInitialize = maxHeight + radius;
xF = (2*radius*rand())-radius;
yF = (2*radius*rand())-radius;
zF = (heightInitialize*rand())-radius;
while(boundaryCondition == true)
n = n + 1;
xI = xF;
yI = yF;
zI = zF;
xF = xI + (2*N*rand())-N;
yF = yI + (2*N*rand())-N;
zF = zI + (2*N*rand())-N;
dist = sqrt(xF^2 + yF^2);
length = sqrt(xF^2 + yF^2 + zF^2);
if (zF >= maxHeight)
flag = 1;
disp(n);
break;
end
if (zF >= 0 && zF < maxHeight && dist < radius)
v1=[xI,yI,zI];
v2=[xF,yF,zF];
v=[v2;v1];
plot3(v(:,1),v(:,2),v(:,3));
pause(pauseTime);
elseif (zF < 0 && length < radius)
v1=[xI,yI,zI];
v2=[xF,yF,zF];
v=[v2;v1];
plot3(v(:,1),v(:,2),v(:,3));
pause(pauseTime);
else
boundaryCondition = false;
while(boundaryCondition == false)
xF = xI + (2*N*rand())-N;
yF = yI + (2*N*rand())-N;
zF = zI + (2*N*rand())-N;
length = sqrt(xF^2 + yF^2 + zF^2);
dist = sqrt(xF^2 + yF^2);
if ((zF >= 0 && zF < maxHeight && dist < radius)||(zF < 0 && length < radius))
boundaryCondition = true;
v1=[xI,yI,zI];
v2=[xF,yF,zF];
v=[v2;v1];
plot3(v(:,1),v(:,2),v(:,3));
pause(pauseTime);
end
end
end
end
anArray(storageCounter)=n;
disp(anArray);
end
plot3(xF,yF,zF,'ok','MarkerFaceColor','r');
hold off
disp(anArray)
disp(n)
function ColorSet=varycolor(NumberOfPlots)
% VARYCOLOR Produces colors with maximum variation on plots with multiple
% lines.
%
% VARYCOLOR(X) returns a matrix of dimension X by 3. The matrix may be
% used in conjunction with the plot command option 'color' to vary the
% color of lines.
%
% Yellow and White colors were not used because of their poor
% translation to presentations.
%
% Example Usage:
% NumberOfPlots=50;
%
% ColorSet=varycolor(NumberOfPlots);
%
% figure
% hold on;
%
% for m=1:NumberOfPlots
% plot(ones(20,1)*m,'Color',ColorSet(m,:))
% end
%Created by Daniel Helmick 8/12/2008
narginchk(1,1)%correct number of input arguements??
nargoutchk(0, 1)%correct number of output arguements??
%Take care of the anomolies
if NumberOfPlots<1
ColorSet=[];
elseif NumberOfPlots==1
ColorSet=[0 1 0];
elseif NumberOfPlots==2
ColorSet=[0 1 0; 0 1 1];
elseif NumberOfPlots==3
ColorSet=[0 1 0; 0 1 1; 0 0 1];
elseif NumberOfPlots==4
ColorSet=[0 1 0; 0 1 1; 0 0 1; 1 0 1];
elseif NumberOfPlots==5
ColorSet=[0 1 0; 0 1 1; 0 0 1; 1 0 1; 1 0 0];
elseif NumberOfPlots==6
ColorSet=[0 1 0; 0 1 1; 0 0 1; 1 0 1; 1 0 0; 0 0 0];
else %default and where this function has an actual advantage
%we have 5 segments to distribute the plots
EachSec=floor(NumberOfPlots/5);
%how many extra lines are there?
ExtraPlots=mod(NumberOfPlots,5);
%initialize our vector
ColorSet=zeros(NumberOfPlots,3);
%This is to deal with the extra plots that don't fit nicely into the
%segments
Adjust=zeros(1,5);
for m=1:ExtraPlots
Adjust(m)=1;
end
SecOne =EachSec+Adjust(1);
SecTwo =EachSec+Adjust(2);
SecThree =EachSec+Adjust(3);
SecFour =EachSec+Adjust(4);
SecFive =EachSec;
for m=1:SecOne
ColorSet(m,:)=[0 1 (m-1)/(SecOne-1)];
end
for m=1:SecTwo
ColorSet(m+SecOne,:)=[0 (SecTwo-m)/(SecTwo) 1];
end
for m=1:SecThree
ColorSet(m+SecOne+SecTwo,:)=[(m)/(SecThree) 0 1];
end
for m=1:SecFour
ColorSet(m+SecOne+SecTwo+SecThree,:)=[1 0 (SecFour-m)/(SecFour)];
end
for m=1:SecFive
ColorSet(m+SecOne+SecTwo+SecThree+SecFour,:)=[(SecFive-m)/(SecFive) 0 0];
end
end
end
