%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Section 5.16.3 Application: Neutralization process
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Plant (5.91): G(s)=-2kd/(1+tauh*s), Gd(s)=kd/(1+tauh*s), kd=2.5e6, tauh=1000
kd=2.5e6;tauh=1000;
G=nd2sys(1,[tauh 1],-2*kd);
Gd=nd2sys(1,[tauh 1],kd);
w=logspace(-4,4,61);
G_f=frsp(G,w);
Gd_f=frsp(Gd,w);

%Figure 5.21
figure(1);clf
wd = [1 2500; 0.05 2500; 2 Inf]
wd2 = [0.02 2500; 0.01 2500; 2 Inf]
vplot('liv,lm',Gd_f,G_f,1,'w:',wd,'w:',wd2,'w:');
xlabel('FREQUENCY [RAD/S]');ylabel('MAGNITUDE');
text(0.1, 3000,'Gd');text(.1, 2e5,'G');
text(1500,0.03,'2500');

% n Buffer tanks
w=logspace(-2,4,121);
tank1=nd2sys(1, [1 1]);
tank1_f = frsp(tank1,w);

%Figure 5.23
figure(2);clf;
vplot('liv,lm',tank1_f);
axis([1.e-1,1e2,1.e-3,1]);
hold;

n=4;
for i=2:n 
tank0=nd2sys(1, [1/i 1]); tank=tank0;
for j=2:i tank=mmult(tank,tank0);end
tank_f = frsp(tank,w);
vplot('liv,lm',tank_f);end
xlabel('Frequency*tauh');ylabel('Magnitude');
text(14,.1,'n=1');text(16,0.02,'n=2');
text(18,0.005,'n=3');text(11,0.002,'n=4');

