matlab程序设计，自己编写的总是不对

在边长L=100的正方形的区域内随机投放半径r=1的圆盘颗粒，保证每个圆盘颗粒不重叠(相交)，圆盘的堆积密度A=0.3,如何程序实现？

n=1;m=1;
C=n*pi;
A(1)=2;B(1)=2;
P(1)=2;Q(1)=2;
while C<3000
m=m+1;
P(m)=rand(1,1)*98+1;Q(m)=rand(1,1)*98+1;
T=1;sum(1)=2
for z=1;n
while T==1
D(z)= A(z)+B(z)*i;
E(m)= P(m)+Q(m)i;
sum(z)=abs(E(m)-D(z));
if sum(z)<2
T=2;
end
end
end
if T==1
n=n+1;
A(m)=P(n);B(m)=Q(n);
C=npi;
end
end

你好同学，我取1000个颗粒循环1000次左右就没有什么重叠了，你可以仔细看看，代码供参考：

function pos = randGenerator()
kn = 1000; xmin = 0; xmax = 100; ymin=0; ymax=100;
r = 1; mass = pi*r^2*1*3000;
np = 1000;% 颗粒个数，自己取
coor = 1+98*rand(np,2);
cn = 10;
pos = coor;
radius = ones(np,1);
velhalf = zeros(size(coor));
velhalfold = velhalf;
dis = velhalf;
dt = 1e-2;
nt = 1200;%循环次数
for tt = 1:1:nt
    fprintf('tt=%d\n',tt);
    f = getForceState(pos,radius, kn, xmin, xmax, ymin, ymax);
    velhalf =  ((2-cn*dt)*velhalfold + 2*dt*f/mass)/(2+cn*dt);
    dis = dis + velhalf*dt;
    pos = coor + dis;
    velhalfold = velhalf;
    if(mod(tt,100)==0)
        figure(1);clf
        showGrains(pos, radius)
        axis([0,100,0,100])
        pause(0.0001)
    end
end
end

function f = getForceState(p, r, kn, xmin, xmax, ymin, ymax)
n = size(p,1);
f1 = zeros(n,2);
for i = 1:1:n-1
    for j = i+1:1:n
        rij = p(j,:)-p(i,:);
        dij = norm(rij);
        dir = rij/dij;
        if(dij<(r(i)+r(j)))
            delta = (r(i)+r(j)) - dij;
            f1(i,:) = f1(i,:) - 100*kn*delta/(r(i)+r(j))*dir;
            f1(j,:) = f1(j,:) + 100*kn*delta/(r(i)+r(j))*dir;
        end
    end
end
f2 = zeros(n,2);
kn = kn*100;
for i = 1:1:n
    if(p(i,1)-xmin< r(i))
        f2(i,1) = f2(i,1) +  kn*(r(i)-(p(i,1)-xmin))/r(i);
    end
    if(-p(i,1)+xmax< r(i))
        f2(i,1) = f2(i,1) -  kn*(r(i)-(-p(i,1)+xmax))/r(i);
    end
    
    if(p(i,2)-ymin< r(i))
        f2(i,2) = f2(i,2) +  kn*(r(i)-(p(i,2)-ymin))/r(i);
    end
    if(-p(i,2)+ymax< r(i))
        f2(i,2) = f2(i,2) -  kn*(r(i)-(-p(i,2)+ymax))/r(i);
    end
end
f=f1+f2;
end

function showGrains(coor, rad)
np = size(coor,1);
col = jet(np);
n = 20;
theta = linspace(0,2*pi,n+1);
x = cos(theta); y = sin(theta);
hold on
for i = 1:1:np
    fill(coor(i,1)+rad(i)*x, coor(i,2)+rad(i)*y, col(i,:));
end
axis('equal')
end

1200次效果（随机的）：

答题不易，如有帮助，还请给个采纳支持一下答主哦