% Collisions script
% A collision between two masses.  1D.
% Andy Ruina  March 17, 2021


% Let Matlab do the algebra for you.
% Look up any commands used here that you don't know.

% Two masses:  A and B.
% Two times:   1)some time  before collision and 
%              2)some time  after  collision.
%Assume initial positions and times  are such that the collision happens.

%All variables are declared here as sybolic and real 
syms mA mB   vA1 vB1  vA2 vB2  e   L1 L2   EK1  EK2 Impulse loss real
 
%Assume you are given mA, mB, vA1, vB1 and e.  Find all else.
L1 = mA*vA1 + mB*vB1;        %total system momentum before collision
L2 = mA*vA2 + mB*vB2;        %total system momentum after  collision

%EVERYTHING INTERESTING IS EQ1 adn EQ2 below. All else is simpler.
% The change in mA momentum is equal and opposite to the mB's change
eq1 =  L1 == L2;             %eq1: is Conservation of momentum
eq1 =  L1-L2==0; eq1 = L1-L2;      %  Equivalent

% Relative velocity reverses and is attenuated by a factor of  0<e<1
eq2 = (vA2-vB2) == -e*(vA1 - vB1); % restitution equation

%SOLVE. A great Matlab command.  
%       Solve a list of equations for a listof variables. 
[vA2,vB2] =  solve([eq1,eq2], [vA2, vB2]);   % Learn it, love it.

% Calculate various motion quantities that result from the collision.
Impulse = mB*(vB2-vB1);   % Integral of F dt in the collision

EK1 = (mA* vA1^2 + mB* vB1^2)/2; %Total kinetic energy before collision
EK2 = (mA* vA2^2 + mB* vB2^2)/2; %Total kinetic energy after  collision
EK1 = simplify(EK1);  EK2 = simplify(EK2); % When in doubt, simplify.

loss = simplify(EK1-EK2); %Loss of kinetic energy. 

%ALL THE REST, BELOW, IS JUST OUTPUT OF THE ABOVE CALCULATIONS.
%Everything you could want to know about a 1D collision.
% Read and understand the resulting output.




fprintf('---------------------------------------------------\n')
fprintf('Formulas for various results of the collision.\n')
fprintf('   vA2     = %s\n'  , char(vA2)    ); %'char' converts symbolic to 
fprintf('   vB2     = %s\n'  , char(vB2)    ); % charactor. 
fprintf('   Impulse = %s\n'  , char(Impulse));
fprintf('   loss    = %s\n', char(loss)   );
fprintf('             (Note: Loss is always >0, so long as |e|<1)\n\n' );



% An example with numbers
mA = 1; mB = 1;  vA1 = 1; vB1 = 0;  e = 1;
fprintf('Example: 2nd mass initially stationary, e = 1 (An ''elastic'' collision)\n')
fprintf('  vA2 = %6.2f ;   '     ,eval(vA2));  %evaluate symbolic expression
fprintf('  vB2 = %6.2f ;   '     ,eval(vB2)); 
fprintf('  Impulse  = %6.2f ;   ',eval(Impulse)); 
fprintf('  loss = %6.2f ;   \n'  ,eval(loss));
fprintf('    (Note the exchange of velocities.)\n\n')
 
% Another  example
mA = 1; mB = 1;  vA1 = 1; vB1 = 0;  e = 0;
fprintf('Example: 2nd mass initially stationary, e = 0 (A ''plastic'' collision)\n')
fprintf('  vA2 = %6.2f ;   '     ,eval(vA2));  %evaluate symbolic expression
fprintf('  vB2 = %6.2f ;   '     ,eval(vB2)); 
fprintf('  Impulse  = %6.2f ;   ',eval(Impulse)); 
fprintf('  loss = %6.2f ;   \n'    ,eval(loss));
fprintf('    (Note the equal velocities.)\n\n')

% Yet another  example
mA = 1; mB = 1;  vA1 = 1; vB1 = -1;  e = 0;
fprintf('Example: opposite initial velocities, e = 0 (A ''plastic'' collision)\n')
fprintf('  vA2 = %6.2f ;   '     ,eval(vA2));  %evaluate symbolic expression
fprintf('  vB2 = %6.2f ;   '     ,eval(vB2)); 
fprintf('  Impulse  = %6.2f ;   ',eval(Impulse)); 
fprintf('  loss = %6.2f ;   \n'    ,eval(loss));
fprintf('    (The end of the universe.)\n')
fprintf('---------------------------------------------------\n\n')




mA = 1; mB = 1;  vA1 = 1; vB1 = -1;  e = -1;
fprintf('Example: opposite initial velocities, e = -1 (A ''???'' collision)\n')
fprintf('  vA2 = %6.2f ;   '     ,eval(vA2));  %evaluate symbolic expression
fprintf('  vB2 = %6.2f ;   '     ,eval(vB2)); 
fprintf('  Impulse  = %6.2f ;   ',eval(Impulse)); 
fprintf('  loss = %6.2f ;   \n'    ,eval(loss));
fprintf('    (The universe is unchanged.)\n')
fprintf('---------------------------------------------------\n\n')


 
