Here is the code:
%%Elastic Collision
clc
clearvars
close all
%-----------------------
N=5; % number of particles
L=10; %length of the box
dt=0.0003;
%X=zeros(2,N);
%V=zeros(2,N);
%R=0*rand(1,N)+0.2; %radius of each particle
%M=2*R; % mass of each particle
theta=0:2*pi/20:2*pi;
nstep=20000; % number of time steps or duration
pit=200; % number of iterations between two plots or speed
%%---------- manually initilizing position adn velocity
%position
X=[1 5 0 3 7;... % x pos %need to keep 2 row and only expand number of columns
3 2 4 5 -8]; % y pos
%velocity
V=[3 5 2 -10 -3;... % x
0 -5 -2 -3 -2]; % y
%Radius of each particule
R=[0.6 0.6 0.6 0.6 0.6] ;
%Mass of each particle
M=[1 1 1 1 1];
%
par1=V(:,1);
par2=V(:,2);
par3=V(:,3);
%auxiliary variables
Ax=zeros(N,N);
Ay=zeros(N,N);
Ar=zeros(N,N);
for i=1:N
for j=i:N
Ar(i,j)=R(1,i)+R(1,j);
end
end
hold on
axis equal
for i=1:N
x=X(1,i)+R(1,i)*cos(theta);
y=X(2,i)+R(1,i)*sin(theta);
plot(x,y,'color','b')
end
plotdomain(L)
grid on
grid minor
drawnow
pause(1)
for T=0:nstep
X=X+dt*V;
%check if particules collided with eachother
for i=1:N
Ax(i,:)=X(1,i)-X(1,:);
Ay(i,:)=X(2,i)-X(2,:);
end
Ax=triu(Ax);
Ay=triu(Ay);
Nrm=(Ax.^2+Ay.^2).^(0.5)-Ar-10^-5; % calculate the distance between each two particles
Nrm=triu(Nrm(1:N-1,2:N));
[row,col]=find(Nrm<0); % find the particles that collided
l=length(row);
%if particles collided calculate the new velocities
if l~=0
col=col+1;
for t=1:l
i=row(t,1);
j=col(t,1);
C1=(2*M(1,j)/(M(1,i)+M(1,j)))*dot(V(:,i)-V(:,j),X(:,i)-X(:,j))/(norm(X(:,i)-X(:,j))^2);
C2=(2*M(1,i)/(M(1,i)+M(1,j)))*dot(V(:,j)-V(:,i),X(:,j)-X(:,i))/(norm(X(:,i)-X(:,j))^2);
V(:,i)=V(:,i)-C1*(X(:,i)-X(:,j));
V(:,j)=V(:,j)-C2*(X(:,j)-X(:,i));
% change color upon collision % add here??
end
end
for i=1:N
%check if particules collided with side walls
if X(1,i)+R(1,i)>=L-10^-5 || X(1,i)-R(1,i)<=-L+10^-5
V(1,i)=-V(1,i);
end
%%check if particules collided with top and bottom walls
if X(2,i)+R(1,i)>=L-10^-5 || X(2,i)-R(1,i)<=-L+10^-5
V(2,i)=-V(2,i);
end
end
%plot
if T==pit*ceil(T/pit)
clf
hold on
axis equal
for i=1:N
x=X(1,i)+R(1,i)*cos(theta);
y=X(2,i)+R(1,i)*sin(theta);
plot(x,y,'color','b')
end
plotdomain(L)
title(T*dt)
grid on
grid minor
drawnow
pause(0.01)
end
end
function plotdomain(L)
plot([-L,-L],[-L,L],'color',[0 0 0])
plot([L,L],[-L,L],'color',[0 0 0])
plot([-L,L],[-L,-L],'color',[0 0 0])
plot([-L,L],[L,L],'color',[0 0 0])
plot(-L-0.1,-L-0.1,'color',[1 1 1])
plot(L+0.1,L+0.1,'color',[1 1 1])
end