%%% file to generate figures opt_cher91_concb & opt_cher91_alphab
%%% arrow.m needed
%%% http://www.mathworks.de/matlabcentral/fileexchange/278-arrow-m
clear; close all;

V0=100; c10=3.3; c20=5.5;
c1f=9.04; c2f=0.64;

%%%%%%%%%%% time optimal diafiltration %%%%%%%%%%%%%%%%%
figure(1)
C1=logspace(log10(10),log10(25),40);
b0=63.42; b1=-12.439; b2=-7.836;
C2=1./exp((b0 + b1 + b2 + b1*log(C1))/b2);
plot(C1,C2,'.','MarkerSize',10); hold on;
axis([0 25 0 6]);

alpha=b1/(b1+b2);
c11=exp(-(b0 + b1 + b2 + b2*log(c20))/b1);
c12=(c11^alpha*c20*c1f*c2f^alpha)/(c20^alpha*c1f^alpha*c2f);
C11=logspace(log10(c11),log10(c12),10);
C22=1./exp((b0 + b1 + b2 + b1*log(C11))/b2);
plot([c10 c11 C11 c1f],[c20 c20 C22 c2f]);

figure(2)
dt1=c10*V0*c20^(b2/b1)/b1*exp(b0/b1)*(mfun('Li',exp(-b0/b1)/c11/c20^(b2/b1))-...
    mfun('Li',exp(-b0/b1)/c10/c20^(b2/b1)));
dt2=c10*V0*(-1/c12+1/c11)/(-b1-b2)/(1-alpha);
h(1)=plot([0;dt1;dt1;dt1+dt2],[0;0;alpha;alpha]); hold on
g=arrow([dt1+dt2;alpha],[dt1+dt2;1.05],'FaceColor',[0 0 1],'EdgeColor',[0 0 1],'length',12,'tipangle',12);

%%%%%%%%% C-CVD %%%%%%
figure(1)
C1=linspace(c10,c1f,10); C2=linspace(c20,c2f,30); 
plot([C1 c1f*ones(size(C2))],[c20*ones(size(C1)) C2],'.-r','MarkerSize',10); hold on;
plot(3.3,5.5,'go','MarkerSize',8);
plot(9.04,0.64,'rx','MarkerSize',10);
xlabel('c1'); ylabel('c2');
legend('S(c_1,c_2)=0','min. time','C-CVD','Location','NorthEast');

dt1=c10*V0*c20^(b2/b1)/b1*exp(b0/b1)*(mfun('Li',exp(-b0/b1)/c1f/c20^(b2/b1))-...
    mfun('Li',exp(-b0/b1)/c10/c20^(b2/b1)));
% dt2=-(V0*(c20^2 - c2f^2))/(2*c2*c1f*(b0 + b2*log(c2) + b1*log(c1f)));

dt2=V0*c10/b2/c1f*(log(b0+b2*log(c20)+b1*log(c1f))-log(b0 + b2*log(c2f) + b1*log(c1f)));

figure(2)
T1=linspace(0,dt1,30); T2=linspace(dt1,dt1+dt2,40);
h(2)=plot([T1 T2],[zeros(size(T1)) ones(size(T2))],'.-r','MarkerSize',10); hold on
axis([0 dt1+dt2+0.1 -.1 1.1]);
xlabel('t'); ylabel('\alpha');
legend(h,'min. time','C-CVD','Location','SouthEast');