Data point in dimension

 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

Answer 1

[Image Analyst]

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