%Simulation of two compartment model
%Data from Riggs p 205 ff creatinin
%   kja 070612

aminit;

k21=0.719;k01=0.461;k12=0.652;k02=0;
A=[-k21-k01 k12;k21 -k21];
C=[1 0];
pd=[-1;-5];
p=[-pd(1)-pd(2);pd(1)*pd(2)];
lam=eig(A);
a=[-lam(1)-lam(2);lam(1)*lam(2)];
Wo=[C;C*A];
Wot=[1 0;a(1) 1];
L=Wo\(Wot\(p-a));
%   Use acker to check results
Lc=acker(A',C',pd)';
%   Form augmented system 
Aa=[A zeros(2,2);L*C A-L*C];
Ca=[C zeros(1,2);zeros(1,2) C];
Ba=zeros(4,1);
sys=ss(Aa,Ba,Ca,0);
%   Simulate
x0=[1;0;0;0];tfinal=6;
[y,t,x]=initial(sys,x0,tfinal);

subplot(221)
pl=plot(t,y(:,1),'b',t,y(:,2),'r--');
set(pl,'Linewidth', AM_data_linewidth);
amaxis([0 6 0 1]);
xlabel('t'); ylabel('x');
lh = legend('actual', 'estimated', 'Location', 'Northeast');
lh = legend(lh, 'boxoff');

subplot(221); hold on;
pl=plot(t,x(:,2),'b',t,x(:,4),'r--');
set(pl,'Linewidth', AM_data_linewidth);
amaxis([0 6 0 0.7]);

amtext(3, 0.12, 'c1');
amtext(4.5, 0.25, 'c2');

amprint('compartment-sim.eps');
