%This function takes in positions, velocities, and masses of 
%objects and checks to see if any are colliding. If so,
%it calculates the new velocities based on coservation
%of momentum along the collision axis.

%Henry Peterson 2016

function [VX,VY] = collisions(XX,YY,VX,VY,MS)

%the buffer between particles to determine collisions
buff = 1e-1;

%number of particles 
nParts = length(XX);

%an array to keep track of whick particles have been used
used = zeros(nParts,1);
%an index for the used array
I = 1;
for i=3:nParts-1        %skip the two planets
    for j=(i+1):nParts
        if find(used == j)
            %already collided with one, move on
            continue;
        end
        %colliding
        if abs(XX(i) - XX(j))<buff && abs(YY(i) - YY(j))<buff
            %add it to used array
            used(I) = j;
            I = I + 1;
            
            %here we want to switch components along collision:
            %the vector along collision
            r = [XX(i)-XX(j), YY(i)-YY(j)];  
            
            %compositions of each particle's velocity vector onto r
            compi = (VX(i)*r(1) + VY(i)*r(2))/norm(r);  
            compj = (VX(j)*r(1) + VY(j)*r(2))/norm(r);
            
            %turn compositions into vectors along r
            proji = compi*[r(1)/norm(r),r(2)/norm(r)];  
            projj = compj*[r(1)/norm(r),r(2)/norm(r)];
                      
            %switch projections in such a way that we conserve momentum
            VX(i) = VX(i) - 2*MS(j)/(MS(j)+MS(i)) * ( proji(1) - projj(1) );
            VY(i) = VY(i) - 2*MS(j)/(MS(j)+MS(i)) * ( proji(2) - projj(2) );
            VX(j) = VX(j) - 2*MS(i)/(MS(j)+MS(i)) * (projj(1) - proji(1) );
            VY(j) = VY(j) - 2*MS(i)/(MS(j)+MS(i)) * (projj(2) - proji(2) );
        end
    end
end