%Astrom Hagglund PID
ks=1;
tt=1;
tn=4.463; tu=1.425;
a=ks*tu/tn;
tau=tu/(tu+tn);
ak1 = 3.8 * exp(-8.4* tau + 7.3 *tau^2);
ak2 = 8.4 * exp(-9.6* tau + 9.8 *tau^2);
titd1 = 5.2 * exp(-2.5* tau -1.4 *tau^2);
titd2 =  3.2 * exp(-1.5* tau -0.93 *tau^2);
tit1 = 0.46 * exp(2.8* tau -2.1 *tau^2);
tit2 =  0.28 * exp(3.8* tau -1.6 *tau^2);
tdtd1 = 0.89 * exp(-0.37* tau -4.1 *tau^2);
tdtd2 =  0.86 * exp(-1.9* tau -0.44 *tau^2);
tdt1 = 0.077 * exp(5.0* tau -4.8 *tau^2);
tdt2 =  0.076 * exp(3.4* tau -1.1 *tau^2);
b1 = 0.40 * exp(0.18* tau + 2.8 *tau^2);
b2 = 0.22 * exp(0.65* tau + 0.051 *tau^2);
 
zr1=ak1/a; zr2 = ak2/a;
ti1 = tit1 * tn; ti2 = tit2 * tn; 
td1 = tdt1 * tn; td2 = tdt2 * tn; 

Zr = zr1; Ti = ti1; Td = td1; b = b1;
P = Zr; I = Zr/Ti; D = Zr*Td;
sim('m6ahs');
vys = res;

b = 1;
sim('m6ahs');
vys2 = res;


t=vys(:,1); t2=vys2(:,1);
w=vys(:,2);
u=vys(:,3); u2=vys2(:,3);
y=vys(:,4); y2=vys2(:,4);
subplot(2,1,1)
plot(t,w,'k:',t,y,'k-',t2,y2,'k--');
v=axis;
%axis([v(1) v(2) v(3) 1.5]);
legend(' w', ' y', ' y(b=1)');
ylabel('w,y');

subplot(2,1,2)
plot(t,u,'k-',t2,u2,'k--');
v=axis;
%axis([v(1) v(2) v(3) 1.5]);
legend(' u', ' u(b=1)');
ylabel('u');
xlabel('t');
set(1, 'PaperUnits', 'centimeters')
set(1, 'PaperPosition', [3 10 11 8.25])
print -deps 'obr6ah.eps'
