%第一列
D=60;
xx1=[];
for i=1:6
theta = linspace(0,2*pi,7);
plot(cos(theta)*D-3*D,sin(theta)*D+(i-3)*sqrt(3)*D,'r-');
hold on
plot(-3*D,(i-3)*sqrt(3)*D,'.','markersize',10)
j= 0;
while j< 50
mu = [-3 (i-3)*sqrt(3)];
SIGMA = [1200 0; 0 1200];
x = mvnrnd(mu,SIGMA,1);
if (abs(x(1)) + abs(x(2))/sqrt(3) ) <= 1*D&& abs(x(2)) <= sqrt(3)*D/2
j = j+1;
hold on
x(1)=x(1)-3*D;
x(2)=x(2)+(i-3)*sqrt(3)*D;
xx1=[xx1;x];
plot(x(1),x(2),'r*');
end
end
end
%第二列
xx2=[];
for i=1:6
theta = linspace(0,2*pi,7);
plot(cos(theta)*D-1.5*D,sin(theta)*D+i*sqrt(3)*D-3.5*sqrt(3)*D,'r-');
hold on
plot(-1.5*D,i*sqrt(3)*D-3.5*sqrt(3)*D,'.','markersize',10)
j= 0;
while j< 50
mu = [-1.5 i*sqrt(3)-3.5*sqrt(3)];
SIGMA = [1200 0; 0 1200];
x = mvnrnd(mu,SIGMA,1);
if (abs(x(1)) + abs(x(2))/sqrt(3) ) <= 1*D&& abs(x(2)) <= sqrt(3)*D/2
j = j+1;
hold on
x(1)=x(1)-1.5*D;
x(2)=x(2)+i*sqrt(3)*D-3.5*sqrt(3)*D;
xx2=[xx2;x];
plot(x(1),x(2),'r*');
end
end
end
%第三列
xx3=[];
for i=1:6
theta = linspace(0,2*pi,7);
plot(cos(theta)*D,sin(theta)*D+(i-3)*sqrt(3)*D,'r-');
hold on
plot(0,(i-3)*sqrt(3)*D,'.','markersize',10)
j= 0;
while j< 50
mu = [0 (i-3)*sqrt(3)];
SIGMA = [1200 0; 0 1200];
x = mvnrnd(mu,SIGMA,1);
if (abs(x(1)) + abs(x(2))/sqrt(3) ) <= 1*D&& abs(x(2)) <= sqrt(3)*D/2
j = j+1;
hold on
x(1)=x(1);
x(2)=x(2)+(i-3)*sqrt(3)*D;
xx3=[xx3;x];
plot(x(1),x(2),'r*');
end
end
end
%第四列
xx4=[];
for i=1:6
theta = linspace(0,2*pi,7);
plot(cos(theta)*D+1.5*D,sin(theta)*D+i*sqrt(3)*D-3.5*sqrt(3)*D,'r-');
hold on
plot(1.5*D,i*sqrt(3)*D-3.5*sqrt(3)*D,'.','markersize',10)
j= 0;
while j< 50
mu = [1.5 i*sqrt(3)-3.5*sqrt(3)];
SIGMA = [1200 0; 0 1200];
x = mvnrnd(mu,SIGMA,1);
if (abs(x(1)) + abs(x(2))/sqrt(3) ) <= 1*D&& abs(x(2)) <= sqrt(3)*D/2
j = j+1;
hold on
x(1)=x(1)+1.5*D;
x(2)=x(2)+i*sqrt(3)*D-3.5*sqrt(3)*D;
xx4=[xx4;x];
plot(x(1),x(2),'r*');
end
end
end
%第五列
xx5=[];
for i=1:6
theta = linspace(0,2*pi,7);
plot(cos(theta)*D+3*D,sin(theta)*D+(i-3)*sqrt(3)*D,'r-');
hold on
plot(3*D,(i-3)*sqrt(3)*D,'.','markersize',10)
j= 0;
while j< 50
mu = [3 (i-3)*sqrt(3)];
SIGMA = [1200 0; 0 1200];
x = mvnrnd(mu,SIGMA,1);
if (abs(x(1)) + abs(x(2))/sqrt(3) ) <= 1*D&& abs(x(2)) <= sqrt(3)*D/2
j = j+1;
hold on
x(1)=x(1)+3*D;
x(2)=x(2)+(i-3)*sqrt(3)*D;
xx5=[xx5;x];
plot(x(1),x(2),'r*');
end
end
end
axis square
%第六列
xx6=[];
for i=1:6
theta = linspace(0,2*pi,7);
plot(cos(theta)*D+4.5*D,sin(theta)*D+i*sqrt(3)*D-3.5*sqrt(3)*D,'r-');
hold on
plot(4.5*D,i*sqrt(3)*D-3.5*sqrt(3)*D,'.','markersize',10)
j= 0;
while j< 50
mu = [4.5 i*sqrt(3)-3.5*sqrt(3)];
SIGMA = [1200 0; 0 1200];
x = mvnrnd(mu,SIGMA,1);
if (abs(x(1)) + abs(x(2))/sqrt(3) ) <= 1*D&& abs(x(2)) <= sqrt(3)*D/2
j = j+1;
hold on
x(1)=x(1)+4.5*D;
x(2)=x(2)+i*sqrt(3)*D-3.5*sqrt(3)*D;
xx6=[xx6;x];
plot(x(1),x(2),'r*');
end
end
end
% %第七列
% xx7=[];
% for i=1:10
% theta = linspace(0,2*pi,7);
% plot(cos(theta)*D+6*D,sin(theta)*D+(i-3)*sqrt(3)*D,'r-');
% hold on
% plot(6*D,(i-3)*sqrt(3)*D,'.','markersize',10)
% j= 0;
%
% while j< 10
% mu = [6 (i-3)*sqrt(3)];
% SIGMA = [30 0; 0 30];
% x = mvnrnd(mu,SIGMA,1);
% if (abs(x(1)) + abs(x(2))/sqrt(3) ) <= 1*D&& abs(x(2)) <= sqrt(3)*D/2
% j = j+1;
% hold on
% x(1)=x(1)+6*D;
% x(2)=x(2)+(i-3)*sqrt(3)*D;
% xx7=[xx7;x];
% plot(x(1),x(2),'r*');
% end
% end
% end
% axis square
% %第八列
% xx8=[];
% for i=1:10
%
% theta = linspace(0,2*pi,7);
% plot(cos(theta)*D+7.5*D,sin(theta)*D+i*sqrt(3)*D-3.5*sqrt(3)*D,'r-');
% hold on
% plot(7.5*D,i*sqrt(3)*D-3.5*sqrt(3)*D,'.','markersize',10)
% j= 0;
% while j< 10
% mu = [7.5 i*sqrt(3)-3.5*sqrt(3)];
% SIGMA = [30 0; 0 30];
% x = mvnrnd(mu,SIGMA,1);
% if (abs(x(1)) + abs(x(2))/sqrt(3) ) <= 1*D&& abs(x(2)) <= sqrt(3)*D/2
% j = j+1;
% hold on
% x(1)=x(1)+7.5*D;
% x(2)=x(2)+i*sqrt(3)*D-3.5*sqrt(3)*D;
% xx8=[xx8;x];
% plot(x(1),x(2),'r*');
% end
% end
% end
% %第九列
% xx9=[];
% for i=1:10
% theta = linspace(0,2*pi,7);
% plot(cos(theta)*D+9*D,sin(theta)*D+(i-3)*sqrt(3)*D,'r-');
% hold on
% plot(9*D,(i-3)*sqrt(3)*D,'.','markersize',10)
% j= 0;
%
% while j< 10
% mu = [9 (i-3)*sqrt(3)];
% SIGMA = [30 0; 0 30];
% x = mvnrnd(mu,SIGMA,1);
% if (abs(x(1)) + abs(x(2))/sqrt(3) ) <= 1*D&& abs(x(2)) <= sqrt(3)*D/2
% j = j+1;
% hold on
% x(1)=x(1)+9*D;
% x(2)=x(2)+(i-3)*sqrt(3)*D;
% xx9=[xx9;x];
% plot(x(1),x(2),'r*');
% end
% end
% end
% axis square
% %第十列
% xx10=[];
% for i=1:10
%
% theta = linspace(0,2*pi,7);
% plot(cos(theta)*D+10.5*D,sin(theta)*D+i*sqrt(3)*D-3.5*sqrt(3)*D,'r-');
% hold on
% plot(10.5*D,i*sqrt(3)*D-3.5*sqrt(3)*D,'.','markersize',10)
% j= 0;
% while j< 10
% mu = [10.5 i*sqrt(3)-3.5*sqrt(3)];
% SIGMA = [30 0; 0 30];
% x = mvnrnd(mu,SIGMA,1);
% if (abs(x(1)) + abs(x(2))/sqrt(3) ) <= 1*D&& abs(x(2)) <= sqrt(3)*D/2
% j = j+1;
% hold on
% x(1)=x(1)+10.5*D;
% x(2)=x(2)+i*sqrt(3)*D-3.5*sqrt(3)*D;
% xx10=[xx10;x];
% plot(x(1),x(2),'r*');
% end
% end
% end