%Time and frequency responses to load disturbances and measurement noise
%   kja 070703

aminit;					% initialize MATALB
%Simulation of the six step responses for the basic loop
%The simulation is based on w
%PID Control of a process with 3 equal lags
P = tf(1,[1 3 3 1]); 
%Controller parameters
k=2.0
ki=1.5;
kd=2.0;
T=0.5;
Tf=0.1;
%---------------------------------------------------
%Prels
denf=[Tf^2/2 Tf 1];
numc=[kd k ki];
denc=conv([1 0],denf)
nump=1;
denp=[1 3 3 1];
dencl=conv(denp,denc)+conv([0 0 0 0 1],conv(nump,numc));
denpr=[T^3 3*T^2 3*T 1];
%---------------------------------------------------
%The six systems
sys11=tf(conv(nump,numc),dencl);
sys12=tf(conv(nump,denc),dencl);
sys21=tf(conv(denp,numc),dencl);
sys22=tf(conv(denp,denc),dencl);
sys11f=tf(1,denpr);
sys21f=tf(denp,denpr);
P=tf(nump,denp);
C=tf(numc,denc);
Ci=tf(1,numc);
P0=1;
%---------------------------------------------------
%Simulation
t=0:0.01:30;
y11=step(sys11,t);
y12=step(sys12,t);
y21=step(sys21,t);
y22=step(sys22,t);
y11f=step(sys11f,t);
y21f=step(sys21f,t);
sp=step(P,t);
sc=step(C,t);
sci=impulse(Ci,t);
%---------------------------------------------------
%Frequency responses
w=logspace(-1,2,200);
m12=squeeze(bode(sys12,w));
m21=squeeze(bode(sys21,w));
mp=squeeze(bode(P,w));
mc=squeeze(bode(C,w));
one=ones(size(w));
%
%Gain crossover
[xxx,iii]=min(abs(mp.*mc-1));
wgc=w(iii);
%---------------------------------------------------
%Plot results for load disturbance response
figure(1); clf; 

subplot(321)
pl=plot(t,y12);			% removed: t,sci,'r:',t,sp,'r--');
set(pl,'Linewidth',AM_data_linewidth);
amaxis([0 20 -0.2 0.5])
ylabel('xd');xlabel('t');

subplot(323)
pl=loglog(w,m12,'b-');		% removed: w,w/ki,'r:',w,mp,'r--');
set(gca, 'XTick', [0.1 1 10], 'XMinorGrid', 'off');
set(gca, 'YTick', [0.01 0.1 1], 'YMinorGrid', 'off'); 
set(pl,'Linewidth',AM_bode_linewidth);
amaxis([0.1 10 0.01 2]); 
xlabel('w');ylabel('gxd');
amprint('distresp-load.eps');

%---------------------------------------------------
%Plot results for noise response
figure(2); clf; 

subplot(321)
pl=plot(t,y21,'b-');		% removed: t,sc,'r--',t,ones(size(t))/P0,'r:');
set(pl,'Linewidth',AM_data_linewidth);
amaxis([0 2 0 20])
ylabel('un');xlabel('t');

subplot(323)
pl=loglog(w,m21,'b-');		% removed: w,mc,'r--',w,one,'r:');
set(gca, 'XTick', [0.1 1 10 100], 'XMinorGrid', 'off');
set(gca, 'YTick', [1 10 100], 'YMinorGrid', 'off'); 
set(pl,'Linewidth',AM_bode_linewidth);
amaxis([0.1 100 0.9 100]); 
xlabel('w');ylabel('gun');
amprint('distresp-noise.eps');

