function[]=Billiards(NumberOfBalls)
% Call: Billiards(Number of Balls), e.g.:
% Billiards(10)
if nargin == 0
NumberOfBalls = 4;%规定球的个数
end
close all;
hold on;
drawflag=1;
factor=82; %Adjust this factor when scaling
DT=5e-3;
Bound=[-4 4 -2 2];
BallColour=[[1 0 0];[1 0 0.5];[1 0.5 0];[0 1 0];[0 0 1];[1 1 0];[1 1 1];...
[0 0.3 0];[0 0 0];[0.65 0.65 0.65];[0 0.75 0.75];[0.3 0 0.6];[0.95 0.65 0.75];...
[0.5 0.25 0];[0 0.2 0.4];[0.9 0.4 0.7];[0.4 0.2 0.3];[0.65 0.55 0.15];[0.25 0.35 0.25];[0.5 0 0]];
TableColour=[.4 .5 .8];
%Plothandle
axis(Bound);
set(gca,'Color',TableColour,'xcolor',TableColour,'ycolor',...
TableColour,'PlotBoxAspectRatio',[1 abs((Bound(3)-Bound(4))/(Bound(2)-Bound(1))) 1],'xtick',[],'ytick',[])
%===============================================================
r=0.15+0.25*rand(NumberOfBalls,1);
%Radii==========================================================
for count1=1:NumberOfBalls;
for count2=1:count1,
rmatrix(count2,count1)=r(count1)+r(count2);
end;
end;
rmatrix=rmatrix.*triu(abs(-1+eye(size(rmatrix))));
%Mass ==========================================================
massa=pi*r.^2;
X=[(Bound(2)-Bound(1)-2*max(r))*rand(NumberOfBalls,1)+Bound(1)+max(r),...
(Bound(4)-Bound(3)-2*max(r))*rand(NumberOfBalls,1)+Bound(3)+max(r)];
for j=1:NumberOfBalls;
for i=1:j;
distmatrix(i,j)=sqrt((X(j,1)-X(i,1))^2+(X(j,2)-X(i,2))^2);
end;
end;
%Initial Edgedetectionmatrix====================================
Botsmatrix=(distmatrix-rmatrix)+tril(abs(-1+eye(size(distmatrix))))+eye(size(distmatrix));
while find(Botsmatrix<=0);
X=[(Bound(2)-Bound(1)-2*max(r))*rand(NumberOfBalls,1)+Bound(1)+max(r),...
(Bound(4)-Bound(3)-2*max(r))*rand(NumberOfBalls,1)+Bound(3)+max(r)];
for j=1:NumberOfBalls;
for i=1:j;
distmatrix(i,j)=sqrt((X(j,1)-X(i,1))^2+(X(j,2)-X(i,2))^2);
end;
end;
Botsmatrix=(distmatrix-rmatrix)+tril(abs(-1+eye(size(distmatrix))))+eye(size(distmatrix));
end
%Initial velocities
V=5*(-1+2*rand(NumberOfBalls,2));
%Plot startingpositions
for k=1:NumberOfBalls;
h(k) = plot(X(k,1),X(k,2),'o', 'MarkerEdgeColor',BallColour(mod(k-1,length(BallColour))+1,:),...
'MarkerFaceColor',BallColour(mod(k-1,length(BallColour))+1,:),'MarkerSize',factor*r(k));
end
pause(1e-1);
%Loop
while drawflag==1;
for k=1:NumberOfBalls;
delete(h(k));
end
%Edgedetecton positive======================================
for j=1:NumberOfBalls;
for i=1:2;
d=X(j,i)+r(j)-Bound(2*i);
if d>=0
dt=d/V(j,i);
X(j,i)=X(j,i)-V(j,i)*dt;
V(j,i)=-V(j,i);
end
end
end
%Edgedetecton negative======================================
for j=1:NumberOfBalls;
for i=1:2;
d=X(j,i)-r(j)-Bound(2*i-1);
if d<=0
dt=d/V(j,i);
X(j,i)=X(j,i)-V(j,i)*dt;
V(j,i)=-V(j,i);
end
end
end
%Distancematrix============================================
for j=1:NumberOfBalls;
for i=1:j;
distmatrix(i,j)=sqrt((X(j,1)-X(i,1))^2+(X(j,2)-X(i,2))^2);
end;
end;
%Collisiondetectionmatrix=========================================================
Botsmatrix=(distmatrix-rmatrix)+tril(abs(-1+eye(size(distmatrix))))+eye(size(distmatrix));
%===========================================================================
if find(Botsmatrix<0);
[I,J]=find(Botsmatrix<0);
for i=1:length(I)
normdist=normr([X(I(i),1)-X(J(i),1) X(I(i),2)-X(J(i),2)]);
vaA=(V(I(i),1)*normdist(1)+V(I(i),2)*normdist(2));
vaB=(V(J(i),1)*normdist(1)+V(J(i),2)*normdist(2));
dt=abs(r(I(i))+r(J(i))-distmatrix(I(i),J(i)))/(abs(vaA)+abs(vaB));
X(I(i):J(i),:)=X(I(i):J(i),:)-V(I(i):J(i),:)*dt;
end
for i=1:length(I)
normdist=normr([X(I(i),1)-X(J(i),1) X(I(i),2)-X(J(i),2)]);
vaA=(V(I(i),1)*normdist(1)+V(I(i),2)*normdist(2));
vbA=(-V(I(i),1)*normdist(2)+V(I(i),2)*normdist(1));
vaB=(V(J(i),1)*normdist(1)+V(J(i),2)*normdist(2));
vbB=(-V(J(i),1)*normdist(2)+V(J(i),2)*normdist(1));
va_A=vaA+2*(vaB-vaA)/(1+massa(I(i))/massa(J(i)));
va_B=vaB+2*(vaA-vaB)/(1+massa(J(i))/massa(I(i)));
V(I(i),1)=va_A*normdist(1)-vbA*normdist(2);
V(I(i),2)=va_A*normdist(2)+vbA*normdist(1);
V(J(i),1)=va_B*normdist(1)-vbB*normdist(2);
V(J(i),2)=va_B*normdist(2)+vbB*normdist(1);
X(I(i):J(i),:)=X(I(i):J(i),:)+V(I(i):J(i),:)*dt;
end
end
%Propagation=================================================
X=X+V*DT;
%Plotting====================================================
for k=1:NumberOfBalls;
h(k) = plot(X(k,1),X(k,2),'o', 'MarkerEdgeColor',BallColour(mod(k-1,length(BallColour))+1,:),...
'MarkerFaceColor',BallColour(mod(k-1,length(BallColour))+1,:),'MarkerSize',factor*r(k));
end
drawnow;
end
