Multiple Spheres in 3D Image rather than plot

Hey guys,

I've got a script that will plot multiple spheres given an inputted void fraction and number of spheres. I was wondering if there was a way to do this using a binary image rather than the plot function? In addition what would be the easiest way to log-normally distributed the sphere radius so they are not all the same size? I want to progress onto utilising simulated annealing on the bed of spheres. I've attached my file below. Thanks.

% Prompts the user to enter the inputs

V = input('Enter the void fraction: ');

n = input('Enter the number of particles: ');

% Calculates the total volume of all the spheres, the volume of the cube

%required based on given void fraction, the axis length based on cube

%volume and sets up a matrix array of 'n' columns full of random numbers

%and defines x, y, z, as starting coordinates for a sphere

volume_spheres = (4/3)*pi()*n;

volume_cube = volume_spheres/V;

axis_length = volume_cube^(1/3);

N = axis_length*rand(3, n);

axis([0 axis_length 0 axis_length 0 axis_length]);

hold on;

title('Particle Modeling');

xlabel('X Direction');

ylabel('Y Direction');

zlabel('Z Direction');

[x, y, z] = sphere();

% Instigates a for loop which will continue 'n' number of times, this then

%utilises a while loop which will check, by use of the pythagorean theorem,

%whether each sphere will not intersect with another. If they do intersect

%the coordinate will be replaced with a new random coordinate, and checked

%again.

for i = 1:n; ii = i;

    while i~=1 && ii>1
        ii = ii-1;
        C1 = N(:,i);
        C2 = N(:,ii);

        x1 = C1(1);
        x2 = C2(1);
        y1 = C1(2);
        y2 = C2(2);
        z1 = C1(3);
        z2 = C2(3);

        dx = x1-x2;
        dy = y1-y2;
        dz = z1-z2;

        dis = sqrt(dx.^2+dy.^2+dz.^2);

        if dis<2;
            N(:,i) = axis_length*rand(3,1);
            ii=i;
        else
            N(:,i) = C1;
            N(:,ii) = C2;
        end
    end

    points = N(:,i);
    xx = points(1);
    yy = points(2);
    zz = points(3);

    surf(x+xx, y+yy, z+zz)

end

%The script will now plot the cumulative distribution of the array in the %x, y, and z direction. It will also compare it to the thereotical %distribution which is linear.

hold off xcord = N(1,:); ycord = N(2,:); zcord = N(3,:); theox = [0 axis_length]; theoy = [0 1];

figure(2) title('Cumulative Probability Distribution') subplot(1,3,1) hold on cdfplot(xcord) plot(theox, theoy, 'm--') legend('Actual', 'Theoretical','Location','NorthWest') title('X-Direction') xlabel('X') ylabel('F(X)') hold off

subplot(1,3,2) hold on cdfplot(ycord) plot(theox, theoy, 'm--') legend('Actual', 'Theoretical','Location','NorthWest') title('Y-Direction') xlabel('Y') ylabel('F(Y)') hold off

subplot(1,3,3) hold on cdfplot(zcord) plot(theox, theoy, 'm--') legend('Actual', 'Theoretical','Location','NorthWest') title('Z-Direction') xlabel('Z') ylabel('F(Z)') hold off