%================================================
%drug_standard.m Standard model
%================================================
clear all
%--------------------------------------------------------------
%CHARACTERISTICS:
%
%--------------------------------------------------------------
%Glossary:
% fine    fineness of change in miligramm for equillibrium calc
%--------------------------------------------------------------
%% Define parameters for simulation:
T=6;
D_L=20;
D_0=20;
k=0.06;
fine=0.003; 
D_num=20;
D_me=40;
D_MS=100;
%--------------------------------------------------------------
%Define a time vector:
tn=linspace(0,T,100);
Tn=linspace(0,20,400);
%--------------------------------------------------------------
%Define matrices:
t=[];
D=[];

%--------------------------------------------------------------
%Define functions:

%--------------------------------------------------------------
%% Program
for n=1:D_num
   Dn=(D_L*exp(-(n-1)*k*T)+D_0*((1-exp(-(n-1)*k*T))/(1-exp(-k*T))))*exp(-k*tn);
   t=[t,tn+(n-1)*T];
   D=[D,Dn];
end

Dmen=D_me*(exp(k*Tn)-1);
DMSn=D_MS*(1-exp(-k*Tn));

amount1=length(find(D>=D_me));
amount2=length(find(D>D_MS));
theff=(amount1-amount2)/(length(D)/(D_num*T));
S=['With this dosing regimen the drug level is for ',num2str(theff),' hours therapeutic.'];
disp(S)

%% Plot: --------------------


% ------------------------------------------

%%
figure(4)
hold on
sp=find(diff(DMSn-Dmen>0));
X=[Tn(1:sp(2)),fliplr(Tn(1:sp(2)))];                %#create continuous x value array for plotting
Y=[Dmen(1:sp(2)),fliplr(DMSn(1:sp(2)))];              %#create y values for out and then back
fill(X,Y,'g');                  %#plot filled area
plot(Tn(1:(sp(2)+20)),Dmen(1:(sp(2)+20)),':r',Tn(1:(sp(2)+20)),DMSn(1:(sp(2)+20)))
annotation('textbox','String',{'UNSAFE'});
annotation('textbox','LineStyle','none','String',{'GOOD'});
annotation('textbox','String',{'INEFFECTIVE'});
annotation('textarrow','String',{'UNSAFE & INEFFECTIVE'});
legend('','minimal effective level','max safe level')
xlabel('time interval T in hours')
ylabel('dose D_0 in mg')
annotation('textbox','String',['drug_standard.m acceptable dosing regimens with k=',num2str(k),'; D_me=',num2str(D_me),'; D_MS=',num2str(D_MS)]);
hold off