% Model: Coupled Two-Mass-Model
% 
%             |--k---|   
% x|--c1--m1--|      |--m2--c3--|x
%             |--c2--|
% 
% states    : x1 @ m1 
%           : x2 @ m2

clear
close all

global m1 m2 c1 c2 c3 k
x1_0 = 0.01;
m1 = 1;
m2 = 0.7;
c1 = 1;
c2 = 2;
c3 = 3;
k = 0.01;

A = [0 1 0 0; -c1/m1-c2/m1 -k/m1 c2/m1 k/m1; 0 0 0 1; c2/m2 k/m2 -c2/m2-c3/m2 -k/m2];

lambda = eig(A);
omega0 = abs(lambda);
f0 = omega0/2/pi

% q = [x1 dx1 x2 dx2];
q0 = [x1_0 0 0 0];
[t,q] = ode45(@two_mass,[0 20],q0);

figure
set(gcf,'position',[70 150 560 420])
subplot(2,1,1)
plot(t,q(:,1),t,q(:,3))
ylabel('position')
legend('x1','x2')
subplot(2,1,2)
plot(t,q(:,2),t,q(:,4))
xlabel('time [s]')
ylabel('speed')
legend('dx1','dx2')

figure
set(gcf,'position',[660 150 560 420])
subplot(2,1,1)
plot(q(:,1),q(:,2))
xlabel('position x1')
ylabel('speed dx1')
subplot(2,1,2)
plot(q(:,3),q(:,4))
xlabel('position x2')
ylabel('speed dx2')