%      SyncSpringMass.M
%   
%  Synchronisation of coupled spring-mass system
%  
% J. Lunze
% 
% 13.9.2018
% Version v. 23.12.2018
% for 2nd edition: 17.2.2021
%
echo off
clear
close all

%% Parameters of the spring-mass system
%
%  distance in meters
%  time in seconds
%  mass in kg
N1=4;
c1=0.5;     % kg/s^2 = N/m
c2=0.5;
c3=0.5;
m1=0.3;   % kg
m2=0.3;
m3=0.3;
m4=0.3;
%
A=[0 1; 0 0];
B=[0; 1];
C=[1 0];
D=0;
nTS=2;
L=[c1/m1   -c1/m1        0         0;
      -c1/m2 (c1+c2)/m2  -c2/m2       0;
       0      -c3/m3   (c2+c3)/m3  -c3/m3;
       0        0        -c3/m4    c3/m4];
nGS=4*nTS;
Abar=kron(eye(N1,N1),A)+kron(-L, B*C);
Cbar=kron(eye(N1,N1),C);
SpringMassSyst=ss(Abar, zeros(nGS,1), Cbar, zeros(N1,1));
%
Tend=40;
dt=0.001;
Time=0:dt:Tend;


%%  Behaviour of the system with identical initial velocities
%
figure(1)
x01=[0; 0.5];
x02=x01;
x03=x01;
x04=x01;
X0=[x01; x02; x03; x04];
[Y,Time]=initial(SpringMassSyst, X0, Time);
subplot(4,1,[1 2])
plot(Time, Y(:,1), 'b');
latextitle('$${\tt SyncSpringMass:}$$ Behaviour of the spring-mass system');%
axis([0 Tend 0 20]);
latexylabel('$$y_i$$ in m');
rotateY;
xleer;
yticks([0 10 20]);
%
subplot(4,1,[3 4])
plot(Time, -L*Y', 'b');
axis([0 Tend -0.4 0.4]);
latexylabel('$$u_i$$ in N$$\;$$');
rotateY,
latexxlabel('$$t$$ in s');
xticks([0 10 20 30 40 50]);
yticks([-0.4 0 0.4]);


%%  Behaviour of the system with identical initial velocities
%
figure(2)
x01=[0.8; 0];
x02=[1; 0];
x03=[1.2; 0];
x04=[1.4; 0];
X0=[x01; x02; x03; x04];
Tend2=40;
dt=0.001;
Time2=[0:dt:Tend2];
Y=initial(SpringMassSyst, X0, Time2);
subplot(4,1,[1 2])
plot(Time2, Y(:,1), 'b', Time2, Y(:,3), 'b--');
latextitle('$${\tt SyncSpringMass:}$$ Behaviour of the spring-mass system');%
axis([0 Tend2 0.7 1.5]);
latexylabel('$$y_1, y_4$$ in m');
rotateY;
xleer;
yticks([0.8 1 1.2 1.4]);
%
subplot(4,1,[3 4])
U=-(L*Y')';
plot(Time2, U(:,1), 'b', Time2, U(:,4), 'b--');
axis([0 Tend2 -0.4 0.4]);
latexylabel('$$u_1, u_4$$ in N$$\;$$');
rotateY,
latexxlabel('$$t$$ in s');
xticks([0 10 20 30 40 50]);
yticks([-0.4 0 0.4]);

%%  Synchronisability analysis of a single mass
%
figure(3)
Mass=ss(A, B, C, D);
rlocprintblue(Mass);
latexxlabel('Re');
latexylabel('Im');
rotateY;
axis([-0.5 0.5 -0.5 0.5])
axis('square');
xticks([-0.4 0 0.4]);
yticks([-0.4 0 0.4]);
plot([-6 6],[0 0],'k:', [0 0],[-5 5],'k:');
latextitle('$${\tt SyncSpringMass:}$$ Root locus of a single mass');


%%   Damped masses
%
figure(4)
d1=0.1;
d2=0.1;
d3=0.1;
d4=0.1;
Ad=[0 1; 0 -d1/m1];
Bd=[0; 1];
Cd=[1 0];
Dd=0;
dampedMass=ss(Ad, Bd, Cd, Dd);
rlocprintblue(dampedMass);
latexxlabel('Re');
latexylabel('Im');
rotateY;
axis([-0.5 0.5 -0.5 0.5])
axis('square');
xticks([-0.4 0 0.4]);
yticks([-0.4 0 0.4]);
plot([-6 6],[0 0],'k:', [0 0],[-5 5],'k:');
latextitle('$${\tt SyncSpringMass:}$$ Root locus of a damped mass');
hold off



%%    Damped spring-mass system
%
figure(5)
Abard=kron(eye(N1,N1),Ad)+kron(-L, Bd*Cd);
Cbard=kron(eye(N1,N1),Cd);
dampedSpringMassSyst=ss(Abard, zeros(nGS,1), Cbard, zeros(N1,1));
x01=[0.8; 0];
x02=[1; 0];
x03=[1.2; 0];
x04=[1.4; 0];
X0=[x01; x02; x03; x04];
[Y,Time]=initial(dampedSpringMassSyst, X0, Time);
subplot(2,1,1)
plot(Time, Y(:,1), 'b', Time, Y(:,4), 'b--');
latextitle('$${\tt SyncSpringMass:}$$ Behaviour of the spring-mass system');%
axis([0 Tend 0.8 1.4]);
latexylabel('$$y_1, y_4$$ in m');
rotateY;
xleer;
yticks([0.8 1 1.2 1.4]);
%
subplot(2,1,2)
plot(Time, -L*Y', 'b');
axis([0 Tend -0.4 0.4]);
latexylabel('$$u_i$$ in N$$\;$$');
rotateY,
latexxlabel('$$t$$ in s');
xticks([0 10 20 30 40 50]);
yticks([-0.4 0 0.4]);


%%   Figures
%
epsfigc15('SyncSpringMass');
