[交流] 简单的MATLAB几何建模问题已有2人参与

function spheres

n = 100;
r = rand( 1, n ) / 10;

for i = 1 : n
    color_value = rand( 1, 3 ); % 颜色
    position = rand( 1, 3 );  % 球心坐标
     [ x, y, z ] = unitsphere;
     
     x =  r( i ) * x + position( 1 );
     y =  r( i ) * y + position( 2 );
     z =  r( i ) * z + position( 3 );
     
     h = surf( x, y, z);
     set( h, 'EdgeColor', 'none', ...
         'FaceColor', color_value  );
     hold on
end
axis equal

V = 4 * pi * r .^3 / 3;
VSUM = sum( V )

function [ x, y, z ] = unitsphere

r = 1; % 单位圆
n = 20;

theta = ( -n : 2 : n ) / n * pi;
phi = ( -n : 2 : n )' / n * pi / 2;

cosphi = cos( phi );
cosphi( 1 ) = 0;
cosphi( n + 1 ) = 0;
sintheta = sin( theta );
sintheta( 1 ) = 0;
sintheta( n + 1 ) = 0;

x = r * cosphi * cos( theta );
y = r * cosphi * sintheta;
z = r * sin( phi ) * ones( 1, n + 1 );

4楼: Originally posted by 碗碗儿2954 at 2016-01-28 09:14:49
好心人您好，我看了您的程序是可行的，但是我想要通过比较第 i 个球和前面生成的 (i-1)的球半径和与球心距离的大小来避免球体之间的重合，在总体积小于比如 v 的情况下，由于position只是任意的1*3数组，并不能存储 ...

function spheres
clear
clc

n = 1000;
v = 0.5; % 总体积
VSUM = 0;
Position = zeros( n, 3 ); % 预设1000个圆球球心坐标
R = zeros( n, 1 ); % 预设1000个圆球半径
i = 1;
m = 1;
myflag = false;

while ( VSUM < v )
   
    m = m +1; % while 循环计数
   
    r = rand( 1 ) / 5; % 球半径
    position = rand( 1, 3 );  % 球心坐标

    for k = 1 : i % 判断任意球心距离与相应球半径和大小
        D = sqrt( ( position(1) - Position( k, 1 ) ) ^ 2 + ( position(2) - Position( k, 2 ) ) ^ 2 + ( position(3) - Position( k, 3 ) ) ^ 2);
        Rr = r + R( k );
        if D < Rr % 如果任意两球半径和大于相应的球心距离
            myflag = true; % 设定继续寻下一个球半径和球心的标记
            break; % 跳出for循环
        end
    end
   
    if myflag == true
        myflag = false;
        continue % 进行下一轮while循环
    end
   
    Vtemp = VSUM + 4 * pi * r .^3 / 3;
    if Vtemp > v % 比较总体积与设定体积大小，这里需要讨论（有可能这次while循环中随机生成的半径较大，可以设定条件）
        break
    end
    VSUM = VSUM + 4 * pi * r .^3 / 3;
    color_value = rand( 1, 3 ); % 颜色
    [ x, y, z ] = unitsphere;
    x =  r * x + position( 1 );
    y =  r * y + position( 2 );
    z =  r * z + position( 3 );
    h = surf( x, y, z);
    set( h, 'EdgeColor', 'none', ...
         'FaceColor', color_value  );
     hold on
         
     Position( i, : ) =  position; % 记录球心坐标
     R( i ) = r; % 记录球半径
     i = i + 1;

    if ( i > n ) || ( m > 10000) % 当球数目达设定上限或while循环达到上限，退出while循环
        break
    end
end
axis equal

i
m
VSUM

function [ x, y, z ] = unitsphere

r = 1; % 单位圆
n = 20;

theta = ( -n : 2 : n ) / n * pi;
phi = ( -n : 2 : n )' / n * pi / 2;

cosphi = cos( phi );
cosphi( 1 ) = 0;
cosphi( n + 1 ) = 0;
sintheta = sin( theta );
sintheta( 1 ) = 0;
sintheta( n + 1 ) = 0;

x = r * cosphi * cos( theta );
y = r * cosphi * sintheta;
z = r * sin( phi ) * ones( 1, n + 1 );