% Generates Fig 5.4
% Computer Controlled Systems (3rd ed)
% Author: B. Wittenmark
% Last edit: 1997-07-02
% Copyright (c) 1996 by K. J. Åström, B. Wittenmark and
% Department of Automatic Control, Lund Institute of
% Technology, Lund, Sweden

close all; clear all;

%newplot;
figure(1)
set(gcf,'PaperUnits','centimeters','PaperPosition',[2 2 13 10])
set(gcf,'Units','centimeters','position',get(gcf,'PaperPosition'))

lw=2;
om=1;
zeta=0.7;
t1=1;
zero=0;
K=1;
h=0.25;
% Sampling the motor
a=exp(-zeta*om*h);
b=om*h*sqrt(1-zeta*zeta);
re=a*cos(b);
d1=-2*re;
d2=a*a;
ds=1+d1+d2;
a2=exp(-h/t1);
a1=-(1+a2);
b1=h-t1*(1-a2);
b2=t1-a2*(t1+h);
bs=b1+b2;

% Solution of the polynomial identity AR+BS=D
n=b2*b2-a1*b1*b2+a2*b1*b1;
w1=d1-a1;
w2=d2-a2;
w3=w1*b2-w2*b1;
w4=w1*(a2*b1-a1*b2)+w2*b2;
if zero<0.5 t0=ds/b1; else t0=ds/bs; end
if zero<0.5 r1= b2/b1; else r1= b2*w3/n; end
if zero<0.5 s0=(d1-a1)/b1; else s0= w4/n; end
if zero<0.5 s1= (d2-a2)/b1; else s1=-a2*w3/n; end
clf;
%
[t,x,y]=sim('block504',[0, 10]);
subplot(2,2,1);
plot(t,y(:,3),'k--',t,y(:,1),'b','Linew',lw);
axis([0,10,0,1.5]);
set(gca,'FontSize',9,'Fontname','NewCenturySchlbk')
set(gca,'xtick',[0 5 10]);
set(gca,'ytick',[0 1]);
ylabel('Output','Fontname','NewCenturySchlbk','Fontsize',9)
text(-2.5,1.5,'(a)','Fontname','NewCenturySchlbk','Fontsize',9)
%title('Fig 5.4','fontname','NewCenturySchlbk','FontSize',9);
drawnow;
%
subplot(2,2,3);
[ts,ys]=stairs(t,y(:,2));
%plot(t,y(:,2));
plot(ts,ys,'b','Linew',lw);
axis([0,10,-1,2.5]);
set(gca,'FontSize',9,'Fontname','NewCenturySchlbk')
set(gca,'xtick',[0 5 10]);
set(gca,'ytick',[0 1 2]);
ylabel('Input','Fontname','NewCenturySchlbk','Fontsize',9)
xlabel('Time','Fontname','NewCenturySchlbk','Fontsize',9)
drawnow;
%
% Sampling the motor
h=1;
a=exp(-zeta*om*h);
b=om*h*sqrt(1-zeta*zeta);
re=a*cos(b);
d1=-2*re;
d2=a*a;
ds=1+d1+d2;
a2=exp(-h/t1);
a1=-(1+a2);
b1=h-t1*(1-a2);
b2=t1-a2*(t1+h);
bs=b1+b2;

% Solution of the polynomial identity AR+BS=D
n=b2*b2-a1*b1*b2+a2*b1*b1;
w1=d1-a1;
w2=d2-a2;
w3=w1*b2-w2*b1;
w4=w1*(a2*b1-a1*b2)+w2*b2;
if zero<0.5 t0=ds/b1; else t0=ds/bs; end
if zero<0.5 r1= b2/b1; else r1= b2*w3/n; end
if zero<0.5 s0=(d1-a1)/b1; else s0= w4/n; end
if zero<0.5 s1= (d2-a2)/b1; else s1=-a2*w3/n; end
%
[t,x,y]=sim('block504',[0, 10]);
subplot(2,2,2);
plot(t,y(:,3),'k--',t,y(:,1),'r','Linew',lw);
axis([0,10,0,1.5]);
set(gca,'FontSize',9,'Fontname','NewCenturySchlbk')
set(gca,'xtick',[0 5 10]);
set(gca,'ytick',[0 1]);
ylabel('Output','Fontname','NewCenturySchlbk','Fontsize',9)
text(-2.5,1.5,'(b)','Fontname','NewCenturySchlbk','Fontsize',9)
drawnow;
%
subplot(2,2,4);
[ts,ys]=stairs(t,y(:,2));
plot(ts,ys,'r','Linew',lw);
axis([0,10,-1,2.5]);
set(gca,'FontSize',9,'Fontname','NewCenturySchlbk')
set(gca,'xtick',[0 5 10]);
set(gca,'ytick',[0 1 2]);
ylabel('Input','Fontname','NewCenturySchlbk','Fontsize',9)
xlabel('Time','Fontname','NewCenturySchlbk','Fontsize',9)
drawnow;
