function twomasses2()
%Two masses connected by a spring.  1D mechanics.
%A constant force is suddenly applied to mass 2 at t=0.
%Lecture example  Sept 7, 2012 -Andy Ruina

%NOTE:  By adding comments and writing things more clearly
%the small bit of code in lecture looks huge and intimicating.
%But if your read slowly through, you will see that the essential
%code is just a few lines.

% REALLY everything is built around the 2 lines labeled
% **READ HERE FIRST** near the bottom.

%To know if your computer is thinking or not, lines like
%this sometimes help.
disp([' '])
disp(['***************************************'])
disp(['Calculation started at: ', datestr(now)])
tic;

%Assign parameters using a matlab "structure"
%  (a name followed by a dot).
% Assume consistent units, say, m,kg, N, s
p.m  = [1 1]'; % the two masses
p.k  = .5;      % the spring constant
p.L0 = 3;      % rest length of spring
p.F0 = 1;      %

% Assign time span and initial conditions for ODE Solution.

tspan = linspace(-1,10,101); 
x0    = [0 p.L0]';   %spring is initial relaxed
v0    = [0  0  ]';   %initial speed is zero
z0    = [x0;  v0];

%SOLVE ODE using ODE23, ODE45, etc (type >>help ODE23)
%Output is the times t you put in in tspan and
%One row of zarray for the values of z at those times
%The differential equation is in the first argument, ie, rhs.
%You can put "options" where the [] is, but to keep it simple
% this example doesn't use them.
%All of the parameters are passed from the ode23 call, to the
%ODE23 function (which is a fancier version of your mid-point method)
%to your differential eqaution which is in "rhs".
[t zarray] = ode23(@rhs,tspan, z0,[],p);


%Easier to read if we define x1 and x2 (unpack them from zarray)
x1 = zarray(:,1);   x2 = zarray(:,2);


%The plotting is easier to understand if,for both masses, we only
%show how much they moved relative to their start positions. 
x1= x1 - x1(1);   x2 = x2-x2(1); %note: x1(1) = x1(start) because t(1) = start) 



plot(t,x1,'r.', t,x2,'b') %use colors, so you know which is which.
  %Always label your plots adequately!
  title('Motion of the masses - (by Andy Ruina)')
  xlabel('time')
  ylabel('x1 (red ... ) and x2 (blue solid)')
  
%To know if your computer is still thinking or not, lines like
%this  help.
disp(['Calculation ended at:   ', datestr(now)])
dur = toc;
disp(['It took  ' num2str(toc), '  seconds. Happy?'])

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%THE DIFFERENTIAL EQUATION IS IN THIS FUNCION
%GIVEN z (the present state), it calculates zdot.
function zdot =rhs(t,z,p)
%Although ODE23 doesn't care, we think of the list of numbers z as
%having different meanings.  For clarity, unpack z into variables with
%sensible names.
x=z(1:2);
v=z(3:4);

%Only turn on the force for t>=0
F=0;
if t>=0, F=p.F0; else F=0;  end


T = p.k*(x(2)-x(1) - p.L0);  %Spring tension

%**READ HERE FIRST**
%THE ODEs, the core of the problem, is just the next two lines.
%UNDERSTAND THESE FIRST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  xdot = v;                       % both x's
  vdot = [ T  (-T+F) ]'./ p.m;    % both v's,   'a = F/m'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%All ODE23 cares about is how fast the list of z variables are changing.
zdot = [xdot;vdot];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
