[xiezhh]

xyr=[];
n=100000;
while n
x=100*rand;
y=100*rand;
r=3+2*rand;
if ~isempty(xyr)
d1=sqrt((xyr(:,1)-x).^2+(xyr(:,2)-y).^2);
d2=xyr(:,3)+r;
if all(d1>=d2)
if x>=r && x<=100-r && y>=r && y<=100-r
xyr=[xyr;x y r];
end
if x<r
d1=sqrt((xyr(:,1)-x-100).^2+(xyr(:,2)-y).^2);
d2=xyr(:,3)+r;
if y>=r && y<=100-r
if all(d1>=d2)
xyr=[xyr;x y r;x+100 y r];
end
elseif y<r
d3=sqrt((xyr(:,1)-x).^2+(xyr(:,2)-y-100).^2);
d4=xyr(:,3)+r;
if all(d1>=d2) && all(d3>=d4)
xyr=[xyr;x y r;x+100 y r;x y+100 r];
end
else
d3=sqrt((xyr(:,1)-x).^2+(xyr(:,2)-y+100).^2);
d4=xyr(:,3)+r;
if all(d1>=d2) && all(d3>=d4)
xyr=[xyr;x y r;x+100 y r;x y-100 r];
end
end
end
if x>100-r
d1=sqrt((xyr(:,1)-x+100).^2+(xyr(:,2)-y).^2);
d2=xyr(:,3)+r;
if y>=r && y<=100-r
if all(d1>=d2)
xyr=[xyr;x y r;x-100 y r];
end
elseif y<r
d3=sqrt((xyr(:,1)-x).^2+(xyr(:,2)-y-100).^2);
d4=xyr(:,3)+r;
if all(d1>=d2) && all(d3>=d4)
xyr=[xyr;x y r;x-100 y r;x y+100 r];
end
else
d3=sqrt((xyr(:,1)-x).^2+(xyr(:,2)-y+100).^2);
d4=xyr(:,3)+r;
if all(d1>=d2) && all(d3>=d4)
xyr=[xyr;x y r;x-100 y r;x y-100 r];
end
end
end
else
n=n-1;
end
else
xyr=[x y r];
if x<r
xyr=[xyr;x+100 y r];
end
if x>100-r
xyr=[xyr;x-100 y r];
end
if y<r
xyr=[xyr;x y+100 r];
end
if y>100-r
xyr=[xyr;x y-100 r];
end
end
end

t=[0:0.1:2*pi 0]';
x=[ones(size(t)) cos(t)]*xyr(:,[1 3])';
y=[ones(size(t)) sin(t)]*xyr(:,[2 3])';
h=plot(x,y);
hold on
plot(xyr(:,1),xyr(:,2),'.')
set(h,'color','r')
axis equal
axis([0 100 0 100])
hold off