function y = atoms()
n = 25; % no. of atoms
dt = 0.1; % time
ma=1;
mb=1;
x0=[0 0 0 0];
l=200;
p = rand(n,2)*0.5*l;
v = rand(n,2)+.5;
r=0.5;
for tt=1:5000
idx = abs(p(:,2)-l)<=1 ; % added this line
p(idx,2) = l ; % added this line
plot(p(:,1),p(:,2),'O', 'MarkerFaceColor', 'b');
axis([0 l 0 l])
drawnow;
for i = 1 : n
if p(i, 1) < 0
v(i, 1) =- v(i, 1);
end
if p(i, 1) > l
v(i, 1) = -v(i, 1);
end
if p(i, 2) < 0
v(i, 2) = v(i, 2);
end
end
for i = 1 : n
for j = i + 1 : n
d = sqrt((p(i, 1) - p(j, 1)) ^ 2 + (p(i, 2) - p(j, 2)) ^ 2);
if d < r
va = sqrt(v(i,1)^2 + v(i,2)^2);
vb = sqrt(v(j,1)^2 + v(j,2)^2);
gama = atan((p(i,2)-p(j,2))/(p(i,1)-p(j,1)));
alpha = atan2(v(i,1),v(i,2)) - gama;
beta = atan2(v(j,1),v(j,2)) - gama;
[x,fval ] = fsolve(@gh,x0);
v(i,1)=x(1)*cos(x(3)+gama);
v(i,2)=x(1)*sin(x(3)+gama);
v(j,1)=x(2)*cos(x(4)+gama);
v(j,2)=x(2)*sin(x(4)+gama);
p(i,1)=p(i,1)+v(i,1)*dt;
p(i,2)=p(i,2)+v(i,2)*dt;
p(j,1)=p(j,1)+v(j,1)*dt;
p(j,2)=p(j,2)+v(j,2)*dt;
p(i,1)=p(i,1)+v(i,1)*dt;
p(i,2)=p(i,2)+v(i,2)*dt;
p(j,1)=p(j,1)+v(j,1)*dt;
p(j,2)=p(j,2)+v(j,2)*dt;
p(i,1)=p(i,1)+v(i,1)*dt;
p(i,2)=p(i,2)+v(i,2)*dt;
p(j,1)=p(j,1)+v(j,1)*dt;
p(j,2)=p(j,2)+v(j,2)*dt;
p(i,1)=p(i,1)+v(i,1)*dt;
p(i,2)=p(i,2)+v(i,2)*dt;
p(j,1)=p(j,1)+v(j,1)*dt;
p(j,2)=p(j,2)+v(j,2)*dt;
p(i,1)=p(i,1)+v(i,1)*dt;
p(i,2)=p(i,2)+v(i,2)*dt;
p(j,1)=p(j,1)+v(j,1)*dt;
p(j,2)=p(j,2)+v(j,2)*dt;
p(i,1)=p(i,1)+v(i,1)*dt;
p(i,2)=p(i,2)+v(i,2)*dt;
p(j,1)=p(j,1)+v(j,1)*dt;
p(j,2)=p(j,2)+v(j,2)*dt;
end
end
end
p=p+v*dt;
end
function y = gh(x)
e = 1;
y = [ma*va*cos(alpha) + mb*vb*cos(beta) - ma*x(1)*cos(x(3)) - mb*x(2)*cos(x(4));
va*sin(alpha) - x(1)*sin(x(3));
vb*sin(beta) - x(2)*sin(x(4));
x(1)*cos(x(3)) - x(2)*cos(x(4)) - e*(va*cos(alpha) - vb*cos(beta))];
end
end
