clc; clear
%%% file to generate figures alb_c1_ana_num
c10=15; c2f=0.1; V0=0.0666;
c1_f=[80;120;240]; c2_0=[98.35;146.3;194.3];
b1=2.877; b2=1.698e-1; b3=1.874e-02; b4=5.708e-04; b5=-2.394e-04; b6=9.334e-05;

case_ex=1;
c20=c2_0(case_ex); c1f=c1_f(case_ex);

% to plot optimal surface
C1=linspace(10,90,50);
C2=-(b4*C1-(b4^2*C1.^2-4*b5*b6*C1.^2+4*b1*b6).^(1/2))/(2*b6);

% first phase, UF c10->c11
c11=-(b4*c20 - (b4^2*c20^2 - 4*b5*b6*c20^2 + 4*b1*b5)^(1/2))/(2*b5);
dt1=-V0*c10*(b5*c10-b5*c11-(b6*c20^2+b3*c20+b1)/c10+...
    (b6*c20^2+b3*c20+b1)/c11+(log(c10)-log(c11))*(b2+b4*c20));

% second phase, singular c11->c1f; c20->c22
c22=-(b4*c1f - (b4^2*c1f^2 - 4*b5*b6*c1f^2 + 4*b1*b6)^(1/2))/(2*b6);
C12=linspace(c11,c1f,20);
C22=-(b4*C12 - (b4^2*C12.^2 - 4*b5*b6*C12.^2 + 4*b1*b6).^(1/2))/(2*b6);
% numerical integral of c10*V0/c1^2/q(c1,c2(c1))/(1-alpha(c1,c2(c1))
dt2=integral(@(c1) c10*V0./c1.^2/(1./(2*b1+b2*c1+b3*(-(b4*c1-...
    (b4^2*c1.^2-4*b5*b6*c1.^2+4*b1*b6).^(1/2))/(2*b6))))./(1-...
    c1.*(0.5*b4*(-(b4*c1-(b4^2*c1.^2-4*b5*b6*c1.^2+...
    4*b1*b6).^(1/2))/(2*b6))+b5*c1)/b1),c11,c1f,'RelTol',1e-6);

% third phase, CVD c22->c2f
dt3=(V0*c10*((b6*c22^2)/2-(b6*c2f^2)/2+c22*(b3+b4*c1f) -...
    c2f*(b3+b4*c1f)+log(c22/c2f)*(b5*c1f^2+b2*c1f+b1)))/c1f;

figure(1)
plot(C1,C2,'.'); hold on
plot([c10 c11 C12 c1f],[c20 c20 C22 c2f]);
plot(c10,c20,'go','MarkerSize',8);plot(c1f,c2f,'rx','MarkerSize',10);
legend('S(c_1,c_2)=0','Location','SouthWest')
axis([10 90 0 110]);

figure(2)
alpha=(b4*C22+2*b5*C12)./(b4*C12+b4*C22+2*b5*C12+2*b6*C22);
alpha=(0.5*b4*C12.*C22+b5*C12.^2)./(b5*C12.^2+b4*C12.*C22+b6*C22.^2);
alpha=(0.5*b4*C12.*C22+b5*C12.^2)/b1;
for i=1:length(C22)
    T(i)=integral(@(c1) c10*V0./c1.^2/(1./(2*b1+b2*c1+b3*(-(b4*c1-...
    (b4^2*c1.^2-4*b5*b6*c1.^2+4*b1*b6).^(1/2))/(2*b6))))./(1-...
    c1.*(0.5*b4*(-(b4*c1-(b4^2*c1.^2-4*b5*b6*c1.^2+...
    4*b1*b6).^(1/2))/(2*b6))+b5*c1)/b1),c11,C12(i),'RelTol',1e-6);
end
plot([0 dt1 dt1+T dt1+T(end) dt1+T(end)+dt3],[0 0 alpha 1 1])
axis([0 dt1+T(end)+dt3 -0.1 1.1]);