MATLAB Coder demo - "Atoms" Simulation

location:"directory of MATLAB"\toolbox\coder\codegendemos\coderdemo_atoms\run_atoms.m

MATLAB version: MATLAB R2011a

% atoms = manage_collision(atoms,i,j)
% Sub-function to manage collision of two atoms indexed i and j.
function atoms = manage_collision(atoms,i,j)
% Never inline this function in the generated code. Makes the
% code more readable when doing code generation.
coder.inline('never');
% Extract x, y, vx, and vy properties from the atoms.
xi = atoms(i).x;
yi = atoms(i).y;
xj = atoms(j).x;
yj = atoms(j).y;
vxi = atoms(i).vx;
vyi = atoms(i).vy;
vxj = atoms(j).vx;
vyj = atoms(j).vy;
% Compute normal of collision point
nx = xi-xj; 
ny = yi-yj;
n2=nx*nx+ny*ny;
n1=sqrt(n2);
nx = nx/n1;
ny = ny/n1;
% Compute velocities along the normal
vDotNi = vxi*nx+vyi*ny;
vDotNj = vxj*nx+vyj*ny;
% Compute the new velocities for the two collided atoms
new_vxi = 2*vDotNj*nx - vxj;
new_vyi = 2*vDotNj*ny - vyj;
new_vxj = 2*vDotNi*nx - vxi;
new_vyj = 2*vDotNi*ny - vyi;
% If result is finite, update the structure array to reflect the collision
if isfinite(new_vxi) && isfinite(new_vyi) && isfinite(new_vxj) && isfinite(new_vyj)
    atoms(i).vx = new_vxi;
    atoms(i).vy = new_vyi;
    atoms(j).vx = new_vxj;
    atoms(j).vy = new_vyj;
end

I think there is a mistake in this function.

% Compute the new velocities for the two collided atoms
new_vxi = 2*vDotNj*nx - vxj;
new_vyi = 2*vDotNj*ny - vyj;
new_vxj = 2*vDotNi*nx - vxi;
new_vyj = 2*vDotNi*ny - vyi;

Suppose the mass of two atoms are the same then

% Compute the new velocities for the two collided atoms
new_vxi = vxj;
new_vyi = vyj;
new_vxj = vxi;
new_vyj = vyi;
The two atoms exchanged their velocities each other.

Answer 1

[Alexander Bottema]

Yes, you're correct. The computation for the new velocities should be:

% Compute impulse along collision normal
vDotN = (vxi-vxj)*nx+(vyi-vyj)*ny;
Ix = vDotN*nx;
Iy = vDotN*ny;
% Compute new velocities for the two collided atoms
new_vxi = vxi - Ix;
new_vyi = vyi - Iy;
new_vxj = vxj + Ix;
new_vyj = vyj + Iy;