How to generate random spheres for a given volume fraction inside a specified volume

Dear All,

I want to use thecode, which was suggested in the question by Jan: https://uk.mathworks.com/matlabcentral/answers/461635-randomly-generated-spheres-in-a-volume-in-matlab, for a similar problem, however, the differences are as follows:

1. I want to generate spheres with radius from R1 to R2 in the specified volume. It can be a cylinder or a cube (in the code, which I have passed below it is cube, but I do not know how to change it to cylinder).
2. The spheres inside the specified volume should have a specified volume fraction, i.e. 0.4 of V.
3. I want to have csv with the center location and the radius of each sphere (it can be obtain from P and R matrixes).

The problem is that I have used the modified code from Jan, however, changing the volume fraction does nothing. It is always a similar number spheres despite a different volume fraction.

Code:

function P = GetRandomSpheres2(Fract, Width, Radius)

% INPUT:

% nWant: Number of spheres

% Width: Dimension of 3d box as [1 x 3] double vector

% Radius: [nWant x 1] vector, radii of spheres

% OUTPUT:

% P: [nWant x 3] matrix, centers

P = zeros(3000, 3);

R2 = (2 * Radius(:)) .^ 2; % Squared to avoid expensive SQRT

iLoop = 1; % Security break to avoid infinite loop

index = 1;

V = 0;

Vmax = Width(1)*Width(2)*Width(3);

Vfract = Fract * Vmax;

nWant = Vfract / (4/3*pi*min(Radius)^3);

while V <= Vfract && iLoop < 1e6

newP = rand(1, 3) .* (Width - 2 * Radius(index)) + Radius(index);

% Auto-expanding, need Matlab >= R2016b. For earlier versions:

% Dist2 = sum(bsxfun(@minus, P(1:nValid, :), newP) .^ 2, 2);

Dist2 = sum((P(1:index-1, :) - newP) .^ 2, 2);

if all(Dist2 > R2(1:index-1) + R2(index))

% Success: The new point does not touch existing sheres:

P(index, :) = newP;

V = V + (4/3*pi*Radius(index)^3);

index = index + 1;

end

iLoop = iLoop + 1;

end

end

Demo file:

Width = [50, 50, 50];

Fract = 0.7;

Radius = [2,4];

%

R = randi(Radius, nWant, 1); % Or however your radii are defined

P = GetRandomSpheres2(Fract, Width, R);

figure

axes('NextPlot', 'add', ...

'XLim', [0, 100], 'YLim', [0, 100], 'ZLim', [0, 100]);

view(3);

[X, Y, Z] = sphere();

for k = 1:nWant

surf(X * R(k) + P(k, 1), Y * R(k) + P(k, 2), Z * R(k) + P(k, 3));

end

Can you help me with this?

Bet regards