% Drug administration example
% RMM, 13 Sep 06

aminit;					% initialize macros

% Parameter settings for the model (need comments)
k0=0.1;
k1=0.1;
k2=0.5;
b0=1.5;

% Compartment model definition
Ada = [-k0-k1 k1; k2 -k2];
Bda = [b0;0];
%Cda = [1 0];
Cda=[0 1];
sysol = ss(Ada, Bda, Cda, 0);
gain=-Cda*(Ada\Bda);

% Determine the desired eigenvalue locations based on desired dynamics
zeta = 0.9;
omega0 = 2.2/(10);
lamr = -2*zeta*omega0/2
lami = sqrt(4*zeta^2*omega0^2 - 4*omega0^2)

% Determine where to set the gains
Kda = place(Ada, Bda, [lamr+lami, lamr-lami]);

% Close the loop
syscl = ss(Ada - Bda*Kda, -Bda / (Cda * ((Ada - Bda * Kda) \ Bda)), Cda, 0);

% Simulation parameters
Tper = 0.1;				% simulation period
Tdur = 1;				% pulse duration
Tmax = 50;				% simulation length

% Set up a pulse train for the open loop system
tvec = 0:Tper:Tmax;
us =0.1*ones(size(tvec));
ps=us;
ps(1:50)=0.30;

% Simulate system for pulse-step and step inputs
[yps, tolv, xps] = lsim(sysol, ps, tvec);
[ys, ts] = step(sysol, Tmax);
ys=0.1*ys;

%Settling time
[xxx,iii]=min(abs(gain*0.1*0.95-ys));
tstep=ts(iii)
[xxx,jjj]=min(abs(gain*0.1*0.95-yps));
tps=tolv(jjj)

figure(1); clf; subplot(431); hold on;
plot(ts, 0.1*gain*ones(size(ts)), 'k-', 'Linewidth', AM_ref_linewidth);
plot(ts, ys, 'b-', 'Linewidth', AM_data_linewidth);
amaxis([0 50 0 2]); box on;
% xlabel('Time {\itt} [min]');
ylabel('Concentration C2');

subplot(434);
plot([0; tolv], [0 us], 'b-', 'LineWidth', AM_data_linewidth);
amaxis([0 50 0 0.4]); box on;
xlabel('Time {\itt} [min]');
ylabel('Input dosage');

amprint('compartment-stepresp-1.eps');

figure(2); clf; subplot(432); hold on;
plot(ts,0.1*gain*ones(size(ts)),'k-', 'Linewidth', AM_ref_linewidth);
plot(tolv, yps,'b-', 'Linewidth', AM_data_linewidth);
amaxis([0 50 0 2]); box on;
% xlabel('Time [min]');
ylabel('Concentration C2');

subplot(435);
plot([0; tolv], [0 ps],'b-', 'LineWidth', AM_data_linewidth);
amaxis([0 50 0 0.4]); box on;
xlabel('Time {\itt} [min]');
ylabel('Input dosage');

amprint('compartment-stepresp-2.eps');
