Data point in dimension

2 次查看（过去 30 天）
Laura
Laura 2014-2-25
编辑： Laura 2014-3-3
tic;
% Generate the same random number for the same seed number.
rand('seed',20001);
% Number of data point
N=40;
% Areal fraction
phi=0.6;
% Dimension of the box
L=14/15; % Length of the box
H=1;% Height of the box
% Diameter of particles
d=sqrt(4*phi*L*H/(pi*N));
% Radius of particles
a=d/2;
for nrun=1:1
nameS=strcat(pdir,'InitCondphi60/','InitCondphi60_',int2str(nrun));
load(nameS);
end
XY0in=[x',y']; % 2 columns matrix
for i= 1:N
while i>1
switch ceil(4*rand)
case 1
epsilon=[ -a/4 0];
case 2
epsilon=[+a/4 0];
case 3
epsilon=[0 -a/4];
case 4
epsilon=[0 +a/4];
end % end switch
% New temporary particle location
Temposition(i,:)=XY0in(i,:)+epsilon;
for j=1:i-1
% Account periodic boundary condition in x-direction
x(j)=Temposition(i,1)-Temposition(j,1);
if abs(x(j))>L/2
x(j)=x(j)-L*sign(x(j));
end
dist(j)=sqrt(x(j).^2+(Temposition(i,2)-Temposition(j,2)).^2);
end
% if the particles are not overlapping, then take it
if min(dist)>=d
XY0in(i,:)=Temposition(i,:) ;
% if the particles are overlapping, then come back to their previous location
elseif min(dist)<d
XY0in(i,:)=XY0in(i,:);
end
end % end while command
end
I have 40 particles in the box and the value of x and y are stored in a matrix of dimension 40 x2 . The first column is x-value and second is y-value.
I want to move the particle(2) by small distance and then calculate the distance between that particle to any particle. If the distance is bigger than or equal to diameter, d. Then I will accept the new position of the particle 2. If the distance is small then d, then I will put it back to the original location. Continue to do that for all 40 particles.
Also, I have to account the periodic boundary condition in the x-direction. Their image will be in the other sides of the box. Note that the box has the dimension from -L/2 to L/2 and the height goes from -H/2 to H/2. It means that the final x and y data will stay in the limit of -L/2 to L/2 for x and -H/2 to H/2 for y.
It seems to me that I did something wrong and the code keeps running without finishing. The temposition location data seems to accumulate each run.
Anybody has any idea what I did wrong according to my explanation above?
Thanks

请先登录，再回答此问题。

采纳的回答

Image Analyst
Image Analyst 2014-2-25
I would think it would go like this, though you didn't attach your data so I can't test it. Look over the code and logic and see if you agree. It's much simpler than what you tried.
tic; % Generate the same random number for the same seed number.
rand('seed',20001);
% Number of data points
N=40;
% Areal fraction
phi=0.6;
% Dimension of the box
L=14/15; % Length of the box
H=1;% Height of the box
% Diameter of particles
d=sqrt(4*phi*L*H/(pi*N));
% Radius of particles
a=d/2;
for nrun=1:1
nameS=strcat(pdir,'InitCondphi60/','InitCondphi60_',int2str(nrun));
load(nameS);
end
XY0in=[x',y']; % 2 columns matrix
for i= 1:N
% For each particle....
switch ceil(4*rand)
case 1
epsilon=[ -a/4 0];
case 2
epsilon=[+a/4 0];
case 3
epsilon=[0 -a/4];
case 4
epsilon=[0 +a/4];
end % end switch
% Get new tentative particle location
Temposition(i,:) = XY0in(i,:) + epsilon;
% Clip to box. Particle must stay in box.
Temposition(i, 1) = max(Temposition(i,1), -L/2);
Temposition(i, 1) = min(Temposition(i,1), +L/2);
Temposition(i, 2) = max(Temposition(i,1), -H/2);
Temposition(i, 2) = min(Temposition(i,1), +H/2);
% Calculate distances to other particles
otherParticles = XY0in; % Initialize
% Remove row for this particular particle we're looking at now.
otherParticles(i,:) = [];
% Now calculate distance to all other particles.
distances = sqrt((Temposition(i,1) - otherParticles(:,1))^2 + (Temposition(i,2) - otherParticles(:,2))^2);
minDistance = min(distances);
if minDistance > d
% It's good, make the change permanent.
XY0in(i,:) = Temposition(i,:);
% Nothing to do if it's less than d,
% just leave XY0in alone!
end
end

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Random Number Generation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by