load o426.dat
t=o426(:,1);
y=o426(:,2);
du = 10;
yn = y/du;
n=length(t);
tn=t(1):1/(10*n):t(end);
yn1=interp1(t,yn,tn,'spline');
yn=yn1;t=tn;
[b,i]=max(diff(yn))
dti=t(i)-t(i-1);
b=b/dti;
ti=t(i);
yi=yn(i);
K=yn(end);
a=yi-b*ti;
tu = -a/b
tn = (K-a)/b - tu
plot(t,yn,'k',[tu,tu+tn],[0,K],'k--',[0 t(end)], [K K], 'k:', [tu+tn, ...
                    tu+tn], [0 K], 'k:')
text(ti,yi,'+', 'HorizontalAlignment', 'center')

ylabel('x_1')
xlabel('t')
set(1, 'PaperUnits', 'centimeters')
set(1, 'PaperPosition', [3 10 11 8.25])
print -deps 'id21exm3a.eps'
pause

nfgstr = [ 1 2  3   4  5  6 ; 0 0.104   0.218   0.319   0.410   0.493; 1  0.368    0.271   0.224   0.195   0.161 ]';
fs = tu/tn;
nstr = sum(nfgstr(:,2)<fs);
fn = nfgstr(nstr,2);
gn = nfgstr(nstr,3);
Tstr = gn * tn;
Dstr = (fs-fn)*tn;

nfgbro = [ 1 2  3   4  5  6; 0.000   0.096   0.192   0.268   0.331 0.385; 1.000   0.500   0.440   0.420   0.410   0.400]';
fs = tu/tn;
nbro = sum(nfgstr(:,2)<fs);
fn = nfgbro(nbro,2);
gn = nfgbro(nbro,3);
Tbro= gn * tn;
Dbro = (fs-fn)*tn;

denstr=1;
for i=1:nstr
  denstr=conv(denstr,[Tstr 1]);
end

gstr=tf(K, denstr)
[yo,to]=step(gstr,t(end));
to = to + Dstr;
i1=find(to<=t(end));
ystr=yo(i1);
tstr=to(i1);

denbro=1;
for i=1:nbro
  denbro=conv(denbro,[Tbro/i 1]);
end

gbro=tf(K, denbro)
[yo,to]=step(gbro,t(end));
to = to + Dbro;
i1=find(to<=t(end));
ybro=yo(i1);
tbro=to(i1);

plot(t,yn,'k',tstr,ystr,'k--',tbro,ybro,'k:')
legend('orig','strejc','broida',4)
ylabel('x_1')
xlabel('t')
set(1, 'PaperUnits', 'centimeters')
set(1, 'PaperPosition', [3 10 11 8.25])
print -deps 'id21exm3b.eps'
