generating random particles in a cylinder

I am trying to generate random particles in a cylinder corresponds to a finite volume fraction. Therefore, each volume fraction results in different number of particles. But, as volume fraction increases, say number of particles increases, the simulation time increases dramatically. For example, for 1100 particles (about volume fraction equal to 0.3) it takes just a minute to generate them in a way that they maintain within the cylinder boundaries and also do not overlap. However, if I want 1400 particles corresponds to volume fraction 0.4 it took about 24 hours and still not finished. Any suggestion is appreciated. The code is as below;

    vol_frac=0.400831;
    lz=0.04;
    lx=lz;ly=lx;
    x(1)=(2*lz-2*R)*rand(1)-lz+R;
    y(1)=(2*sqrt(lz^2-x(1)^2)-2*R)*rand(1)-sqrt(lz^2-x(1)^2)+R;
    z(1)=(2*lz-2*R)*rand(1)-lz+R;
    %no. of particles
    Npart=round(6*0.08^3*vol_frac/4/0.006^3);
    %Npart=1100;
    for i=2:Npart;
    x(i)=(2*lz-2*R)*rand(1)-lz+R;
    y(i)=(2*sqrt(lz^2-x(i)^2)-2*R)*rand(1)-sqrt(lz^2-x(i)^2)+R;
    z(i)=(2*lz-2*R)*rand(1)-lz+R;
     k=i-1;
    while k>=1;
         c=((y(k)-y(i))^2+(x(k)-x(i))^2+(z(k)-z(i))^2)^0.5;
             if c<=2*R 
                 x(i)=(2*lz-2*R)*rand(1)-lz+R;
                 y(i)=(2*sqrt(lz^2-x(i)^2)-2*R)*rand(1)-sqrt(lz^2-x(i)^2)+R;
                 z(i)=(2*lz-2*R)*rand(1)-lz+R;
                 k=i-1;
             else
                 k=k-1;
             end
    end

    end

[Image Analyst]

What is R, and is it a constant? Please add comments. Also you should know that typing control-i in the MATLAB editor will fix up the indenting so that when you paste it in here, it does not look all messed up like it does now.

[hamed]

R is particle radius and equal to 0.003

[hamed]

Actually the number of particles has been limited to 1277. I assume there should be sth wrong with the function rand and the way it generates random numbers. and the figure below shows the distribution of generated x values within the boundaries. As can be seen, the distribution of the random numbers is normal that means the particles are accumulated in the centre of cylinder.

[Guillaume]

The issue is not with rand which does generate uniformly distributed numbers but most likely with your mapping from a cube (which is what is generated by rand) to a cylinder.

Can you comment your code? Explain each equation and particularly the purpose of the while loop. Also which direction is the axis of the cylinder?

I would hope that x, y and z have been predeclared with zero. Otherwise, the resizing / copying with each iteration of the for loop would be a reason for slowdowns.

[hamed]

Thanks for your comments. The axial direction of the cylinder is along the z-axis (lz). Polar coordinate has been assumed in x-y plane. According to generate the position of the particles within a circle in x-y plane, I firstly guess an x value then from this value I calculate the y value in a way with which fall within the boundary of the circle. A margin of R equal to radius of particles exclude from the min and max of random values for x,y and z to avoid overlapping with the boundaries. "While" loop is for checking the overlapping of the particles. In each iteration, "for" loop generate a random position (x, y and z) and the "while" loop check whether this value is overlapping with all previous positions. If so, it would be replaced with another x,y and z and again checked. Unfortunately, I forgot to preallocate a zero values for the x,y and z. However, I assume that it might not the cause root of main problem. Thanks again and if you need any further information just don't hesitate to ask.

[hamed]

I'm still not getting appropriate answer. I think the problem is due to random number rejection by overlapping check which results in predictable random numbers. However, I have revised the previous code due to mapping like this and becomes faster;

%Initial random position of particles (Normal distribution)

    vol_frac=0.400831;
    %no. of particles
    rng shuffle
    % Npart=1100;
      R=0.003;
      Npart=round(6*0.08^3*vol_frac/4/0.006^3);
      x=zeros(1,Npart);y=zeros(1,Npart);z=zeros(1,Npart);
      lz=0.04;
      lx=lz;ly=lx;
      teta=2*pi*rand;
      ri=(lz-R)*sqrt(rand);
      x(1)=ri*cos(teta);
      y(1)=ri*sin(teta);
      z(1)=(2*lz-2*R)*rand-lz+R;
      for i=2:Npart;
          teta=2*pi*rand;
          ri=(lz-R)*sqrt(rand);
          x(i)=ri*cos(teta);
          y(i)=ri*sin(teta);
          z(i)=(2*lz-2*R)*rand-lz+R;
          k=i-1;
      while k>=1;
          c=((y(k)-y(i))^2+(x(k)-x(i))^2+(z(k)-z(i))^2)^0.5;
               if c<2*R
                   teta=2*pi*rand;
                   ri=(lz-R)*sqrt(rand);
                   x(i)=ri*cos(teta);
                   y(i)=ri*sin(teta);
                   z(i)=(2*lz-2*R)*rand-lz+R;
                   k=i-1;
               else
                   k=k-1;
               end
      end

      end

[Denisse Campos Muñoz]

I used your code to generating particles in cylindrical coordinates and that these follow uniform distribution using random ('unif', ..., ...) and would like to know how to graph the particles.

How I can plot the distributed along the center of each particle?