clear; close all
%
%  Demo second law of thermo with "molecules" bouncing around
%
box = 100;     %% size of box edge
xmol = [];
ymol = [];
n = 5;
sp = 4.5;
for i=1:n
    for j=1:n
        xmol = [xmol; box/2-sp*i];
        ymol = [ymol; box/2-sp*j];
    end
end
diam = 3;     %% diameter of molecule
rad = diam/2;
nm = length(xmol);
for i=1:nm
    hm(i) = plot(xmol(i),ymol(i),'o','MarkerSize',7);
    hold on
end
diam = 3;  %% estimated
diam2 = diam^2;
axis equal
axis ( (box/2)*[-1 1 -1 1]);
set(gca,'XTick',[],'YTick',[])
angles = 360*rand(1,nm);
theta = angles*pi/180;
v = 1;
vx = v*cos(theta);
vy = v*sin(theta);
t = 0;
forever = 1;
briefly = 0.005;
dt = 0.1;
disp('Hit ENTER to start: '); pause
buffer = (box-diam)/2;
while forever
    for i=1:nm
        xmol(i) = xmol(i) + dt*vx(i);
        ymol(i) = ymol(i) + dt*vy(i);
        set(hm(i),'XData',xmol,'YData',ymol);
    end
    pause(briefly);
    %
    %  Check for wall collissions
    %
    hitx = find(abs(xmol) > buffer);
    for i=hitx
        xmol(i) = sign(xmol(i))*buffer;
        vx(i) = -vx(i);
    end
    hity = find(abs(ymol) > buffer);
    for i=hity
        ymol(i) = sign(ymol(i))*buffer;
        vy(i) = -vy(i);
    end
    %
    %  Check for collissions betw molecules
    %
    for i=1:nm-1
        for j=i+1:nm
            dist2 = (xmol(i)-xmol(j))^2 + (ymol(i)-ymol(j))^2;
            if dist2 < diam2
                theta_coll = atan2(ymol(i)-ymol(j),xmol(i)-xmol(j));
                theta_i = atan2(vy(i),vx(i));
                theta_j = atan2(vy(j),vx(j));
                theta_i = 2*theta_coll - theta_i;
                vi = sqrt(vx(i)^2 + vy(i)^2);
                vx(i) = v*cos(theta_i);
                vy(i) = v*sin(theta_i);
                vj = sqrt(vx(j)^2 + vy(j)^2);
                vx(j) = v*cos(theta_i);
                vy(j) = v*sin(theta_j);
            end
        end
    end         
end