function [] = lautim
%%% file to generate figures lau_tim
%%% arrow.m needed
%%% http://www.mathworks.de/matlabcentral/fileexchange/278-arrow-m
clc; clear
C1=0:30:1000;
b1=0.9800; b2=0.7647; b3=1.8080e-5;
b4=5.8607e-7; b5=2.5066e-10; b6=4.0600e-11;
C2=(2*b3*b4-2*b6+2*b2*b6+(4*b2^2*b6^2-4*b2*b3*b4*b6+12*b2*b3*b5*b6*C1-...
    8*b2*b6^2+4*b3^2*b4^2-12*b3^2*b4*b5*C1+9*b3^2*b5^2*C1.^2+4*b3*b4*b6-...
    12*b3*b5*b6*C1+4*b6^2).^(1/2)-3*b3*b5*C1)/(6*b3*b6);
figure(1)
plot(C1,C2,'.','MarkerSize',10); hold on

c10=100; c20=4000;
c1f=400; c2f=1000;
x0=[c10;c20;1;0]; opts=odeset('RelTol',1e-6);

alpha_max=[1e5;5;2;1];

%alpha=inf
sim=1; xopt=[1.6594e6;3.5419e6;1];
xL=1e-6*ones(3,1); xU=1e30*ones(3,1);
fminconoptions = optimset('Algorithm','active-set', 'Display', 'none');
xopt=fmincon(@objfun,xopt,[],[],[],[],xL,xU,@confun,fminconoptions);
    
x0sim=x0; Y=[];
for i=1:3
    [t,y]=ode45(@diafilt,[0 1],x0sim,opts,[i;xopt(i)]);
    x0sim=y(end,:); Y=[Y;y];
    if i==2
        C1=y(:,1); C2=y(:,2); T=t*xopt(2)+xopt(1);
    end
end
plot(Y(:,1),Y(:,2));

figure(2)
alpha=(C2*(2*b2*b6-2*b2^2*b6-2*b2*b3*b4)+C2.^2*(2*b3^2*b4-2*b3*b6+8*b2*b3*b6)...
    -6*b3^2*b6*C2.^3+C1*(2*b1*b5-2*b1*b2*b5)+C1.*C2*(3*b1*b3*b5+3*b2*b3*b5)...
    -3*b3^2*b5*C1.*C2.^2)/2./(b5*C1+b6*C2-b2*b5*C1-b3*b4*C2-b2*b6*C2+3*b3*b6*C2.^2+3*b3*b5*C1.*C2);
plot([0;xopt(1);T]/sum(xopt),[0;0;alpha]); hold on;
arrow([1;alpha(end)],[1;alpha(end)+0.1],'Length',7,'Width',0.1,'TipAngle',30,'EdgeColor','b','FaceColor','b');
axis([0 1.05 -0.05 1.05]); xlabel('t'); ylabel('a');

% alpha=5
sim=2;
% xopt=[1.6594e6;3.5419e6;1e2];
xL=1e-3*ones(3,1); xU=1e30*ones(3,1);
xopt=fmincon(@objfun,xopt,[],[],[],[],xL,xU,@confun,fminconoptions);

x0sim=x0;
for i=1:3
    [t,y]=ode45(@diafilt,[0 1],x0sim,opts,[i;xopt(i)]);
    x0sim=y(end,:);
end
figure(1)
plot(y(:,1),y(:,2),':k','LineWidth',2);

%alpha=2
sim=3;
xL=1e-3*ones(3,1); xU=1e30*ones(3,1);
xopt=fmincon(@objfun,xopt,[],[],[],[],xL,xU,@confun,fminconoptions);
    
x0sim=x0;
for i=1:3
    [t,y]=ode45(@diafilt,[0 1],x0sim,opts,[i;xopt(i)]);
    x0sim=y(end,:);
end
plot(y(:,1),y(:,2),'-.g');

%alpha=1
sim=4;
xL=1e-3*ones(3,1); xU=1e30*ones(3,1);
xopt=fmincon(@objfun,xopt,[],[],[],[],xL,xU,@confun,fminconoptions);
    
x0sim=x0;
for i=1:3
    [t,y]=ode45(@diafilt,[0 1],x0sim,opts,[i;xopt(i)]);
    x0sim=y(end,:); Y=[Y;y];
end
plot(y(:,1),y(:,2),'--r')

plot(c10,c20,'go','MarkerSize',8); plot(c1f,c2f,'rx','MarkerSize',10);
axis([0 1000 800 9000])
xlabel('c1'); ylabel('c2');
h=legend('S_{time}(c_1,c_2)=0','min. time','\alpha_{max}=5','\alpha_{max}=2','\alpha_{max}=1');

    function f = objfun(x)
        f=sum(x);
    end

    function [c,ceq] = confun(x)
        x0sim=x0;
        for i=1:3
            [~,y]=ode45(@diafilt,[0 1],x0sim,opts,[i;x(i)]);
            x0sim=y(end,:);
        end
        
        c=[]; ceq=[y(end,1)-400;y(end,2)-1000];
    end

    function f = diafilt(t,x,u)
        c1=x(1); c2=x(2); V=x(3); tf=u(2);
        b1=1; b2=0.7647; b3=1.8080e-5;
        b4=5.8607e-7; b5=2.5066e-10; b6=4.0600e-11;

        R1=b1; R2=b2-b3*c2; q=b4-b5*c1-b6*c2;

        switch u(1)
            case 1
                alpha=0;
            case 2
                alpha=(c2*(2*b2*b6-2*b2^2*b6-2*b2*b3*b4)+...
                    c2^2*(2*b3^2*b4-2*b3*b6+8*b2*b3*b6)-6*b3^2*b6*c2^3+...
                    c1*(2*b1*b5-2*b1*b2*b5)+c1*c2*(3*b1*b3*b5+3*b2*b3*b5)-...
                    3*b3^2*b5*c1*c2^2)/2/(b5*c1+b6*c2-b2*b5*c1-...
                    b3*b4*c2-b2*b6*c2+3*b3*b6*c2^2+3*b3*b5*c1*c2);
            case 3
                alpha=alpha_max(sim);
        end

        f=[c1*q/V*(R1-alpha);
            c2*q/V*(R2-alpha);
            (alpha-1)*q;
            alpha*q]*tf;
    end
end