%  Cervical Syndrome (first part)
%  ------------------------------
%
echo on
diary mmps_hw2a.dry
%
%  Enter constants first
%
k1 = 0.3;
k2 = 8;
k3 = 3;
B1 = 0.8;
B2 = 10;
B3 = 12;
m1 = 1.2;
m2 = 14;
m3 = 3.2;
m4 = 24;
kk = k1 + k2 + k3;
BB = B1 + B2 + B3;
%
%  Compute system matrices
%
A = [  0      0      0      0      1      0      0      0
       0      0      0      0      0      1      0      0
       0      0      0      0      0      0      1      0
       0      0      0      0      0      0      0      1
      -k1/m1  k1/m1  0      0     -B1/m1  B1/m1  0      0
       k1/m2 -kk/m2  k2/m2  k3/m2  B1/m2 -BB/m2  B2/m2  B3/m2
       0      k2/m3 -k2/m3  0      0      B2/m3 -B2/m3  0
       0      k3/m4  0     -k3/m4  0      B3/m4  0     -B3/m4 ];
b = [ 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 1/m4 ];
c = [ 1 -1 0 0 0 0 0 0 ];
d = 0;
%
%  Check eigenvalues
%
eig(A)
%
%  Make a system
%
S = ss(A,b,c,d)
%
%  Simulate the system over 15 time periods
%
t = 0:0.15:15;
u = sin(3*pi*t);
x0 = zeros(8,1);
y = lsim(S,u,t,x0);
plot(t,y)
grid on
title('Cervical Syndrome')
xlabel('Time (sec)')
ylabel('Amplitude')
print -djpeg mmps_hw2a1.jpg
clf reset
%
%  Compute the Bode diagram
%
w = logspace(-2,2);
[mag,phase] = bode(S,w);
magdb = 20*log10(mag);
subplot(2,1,1)
semilogx(w,magdb(:))
grid on
title('Cervical Syndrome')
xlabel('Frequency')
ylabel('Magnitude')
subplot(2,1,2)
semilogx(w,phase(:))
grid on
xlabel('Frequency')
ylabel('Phase')
print -djpeg mmps_hw2a2.jpg
clf reset
%
%  Close the diary and finish up
%
diary off
return
