%  Cervical Syndrome (Second Part)
%  -------------------------------
%
echo on
diary mmps_hw2b.dry
%
%  Enter constants first
%
k1nom = 0.3;
k2 = 8;
k3 = 3;
B1nom = 0.8;
B2 = 10;
B3 = 12;
m1 = 1.2;
m2 = 14;
m3 = 3.2;
m4 = 24;
%
%  First case: k1 and B1 nominal - 50%
%
k1 = 0.5*k1nom;
B1 = 0.5*B1nom;
kk = k1 + k2 + k3;
BB = B1 + B2 + B3;
%
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;
%
%  Make a system
%
S = ss(A,b,c,d);
%
w = logspace(-2,2);
[mag,phase] = bode(S,w);
magdb = 20*log10(mag);
mm = magdb(:);
pp = phase(:);
%
%  Second case: k1 = k1nom - 50%, B1 = B1nom + 50%
%
k1 = 0.5*k1nom;
B1 = 1.5*B1nom;
kk = k1 + k2 + k3;
BB = B1 + B2 + B3;
%
A(5:6,:) = [ -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 ];
S = ss(A,b,c,d);
%
[mag,phase] = bode(S,w);
magdb = 20*log10(mag);
mm = [ mm , magdb(:) ];
pp = [ pp , phase(:) ];
%
%  Third case: k1 = k1nom + 50%, B1 = B1nom - 50%
%
k1 = 1.5*k1nom;
B1 = 0.5*B1nom;
kk = k1 + k2 + k3;
BB = B1 + B2 + B3;
%
A(5:6,:) = [ -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 ];
S = ss(A,b,c,d);
%
[mag,phase] = bode(S,w);
magdb = 20*log10(mag);
mm = [ mm , magdb(:) ];
pp = [ pp , phase(:) ];
%
%  Fourth case: k1 = k1nom + 50%, B1 = B1nom + 50%
%
k1 = 1.5*k1nom;
B1 = 1.5*B1nom;
kk = k1 + k2 + k3;
BB = B1 + B2 + B3;
%
A(5:6,:) = [ -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 ];
S = ss(A,b,c,d);
%
[mag,phase] = bode(S,w);
magdb = 20*log10(mag);
mm = [ mm , magdb(:) ];
pp = [ pp , phase(:) ];
%
%  Compute the minima and maxima
%
mmin = min(mm')';
mmax = max(mm')';
pmin = min(pp')';
pmax = max(pp')';
%
%  Plot the minima and maxima together
%
magdb = [ mmin , mmax ];
phase = [ pmin , pmax ];
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_hw2b.jpg
clf reset
%
%  Close the diary and finish up
%
diary off
return
