function xu
%% Radiopaque/Ethylen-Glycol Separation : 

% % initial and final conditions :
% x0 = [120; 5];     param.c10 = x0(1);  param.c20 = x0(2);   % [g/L]
% xf = [400; 0.1];   param.c1f = xf(1);  param.c2f = xf(2);   % [g/L]
% param.V0 = 4;                                               % [L]

% initial and final conditions :
x0 = [12; 0.5];    param.c10 = x0(1);  param.c20 = x0(2);   % [g/dL]
xf = [40; 0.01];   param.c1f = xf(1);  param.c2f = xf(2);   % [g/dL]
param.V0 = 40;                                              % [dL]

% parameters :
param.K    = 1/(param.c10*param.V0);
param.N    = 100000;
param.b    = [-29.19; 118.1; 0.0073; 0.813];
param.clim = 57.1630;
param.copt = 24.5001;
param.eps  = 1e-3;

% weigth coeficients : 
param.wT = 2;
% param.wD = 1e-3;
param.wD  = 10;

% control :
alpha0   = 0;
alpha1   = 1;
alphainf = 1000e1;

% Tolerances :
AbsTol = 1e-13;
RelTol = 1e-12;

% simulation step : 
krok = 0.5;

%% minimum time : 
x01t      = [param.c10;param.c20;0];
options1t = odeset('Events',@event1t,'AbsTol',AbsTol,'RelTol',RelTol);
tspan1t   = [0:krok:param.N];
[time1t,state1t,te1t] = ode45(@model,tspan1t,x01t,options1t,alpha0,param);

options2t = odeset('Events',@event2t,'AbsTol',AbsTol,'RelTol',RelTol);
x02t      = state1t(end,:);
tspan2t   = [time1t(end):krok:param.N];
[time2t,state2t,te2t] = ode45(@model,tspan2t,x02t,options2t,alpha1,param);

options3t = odeset('Events',@event3t,'AbsTol',AbsTol,'RelTol',RelTol);
x03t      = state2t(end,:);
tspan3t   = [time2t(end):krok:param.N];
[time3t,state3t,te3t] = ode45(@model,tspan3t,x03t,options3t,alpha0,param);

%% minimum diluant : 

x01d      = [param.c10;param.c20;0];
options1d = odeset('Events',@event1d,'AbsTol',AbsTol,'RelTol',RelTol);
tspan1d   = [0:krok:param.N];
[time1d,state1d,te1d] = ode45(@model,tspan1d,x01d,options1d,alpha0,param);           
 
options2d = odeset('Events',@event2d,'AbsTol',AbsTol,'RelTol',RelTol);
x02d      = state1d(end,:);
tspan2d   = [time1d(end):krok:param.N];
[time2d,state2d,te2d] = ode45(@model,tspan2d,x02d,options2d,alpha1,param);

options3d = odeset('Events',@event3d,'AbsTol',AbsTol,'RelTol',RelTol);
x03d      = state2d(end,:);
tspan3d   = [time2d(end):krok:param.N];
[time3d,state3d,te3d] = ode45(@model,tspan3d,x03d,options3d,alphainf,param);

%% multi-objective : 

x01m      = [param.c10;param.c20;0];
options1m = odeset('Events',@event1m,'AbsTol',AbsTol,'RelTol',RelTol);
tspan1m   = [0:krok:param.N];
[time1m,state1m,te1m] = ode45(@model,tspan1m,x01m,options1m,alpha0,param);           

options2m = odeset('Events',@event2m,'AbsTol',AbsTol,'RelTol',RelTol);
x02m      = state1m(end,:);
tspan2m   = [time1m(end):krok:param.N];
[time2m,state2m,te2m] = ode45(@model,tspan2m,x02m,options2m,alpha1,param);

options3m = odeset('Events',@event3m,'AbsTol',AbsTol,'RelTol',RelTol);
x03m      = state2m(end,:);
tspan3m   = [time2m(end):krok:param.N];
[time3m,state3m,te3m] = ode45(@model,tspan3m,x03m,options3m,alphainf,param);

%% Figures :
line        = 1.5;
marker_size = 9;
figure
plot([state1d(:,1);state2d(:,1);state3d(:,1)],...
     [state1d(:,2);state2d(:,2);state3d(:,2)],'g--','linewidth',line);       % minimum diluant
hold on      
plot([state1m(:,1);state2m(:,1);state3m(:,1)],...
     [state1m(:,2);state2m(:,2);state3m(:,2)],'r','linewidth',line);       % multi-objective 
plot([state1t(:,1);state2t(:,1);state3t(:,1)],...
     [state1t(:,2);state2t(:,2);state3t(:,2)],'b-.','linewidth',line);       % minimum time
plot(param.c10,param.c20,'go',param.c1f,param.c2f,'rx',...
     'markersize',marker_size)
xlabel('c1')
ylabel('c2')
axis([10 60 -0.05 0.66])
legend('minimum diluant','economically optimal','minimum time',...
       'Location','NorthEast')
% fs=6; fz=11;
% exportfig(gcf,'RaEtCl_state.eps','width',fs,'fontmode','fixed',...
%          'fontsize',fz,'Color','rgb'); 
set(gcf, 'PaperUnits', 'centimeters') 
set(gcf, 'PaperPosition', [3 10 16 10])
print -depsc 'RaEtCl_state.eps'
     
figure
a0d     = alpha0*ones(length(time1d),1);
a1d     = alpha1*ones(length(time2d(1:25)),1);
plot([time1d;time2d(1:25)],[a0d;a1d],'g--','linewidth',line)
hold on
a0m     = alpha0*ones(length(time1m),1);
a1m     = alpha1*ones(length(time2m),1);
plot([time1m;time2m],[a0m;a1m],'r','linewidth',line)
hold on
a0t     = alpha0*ones(length(time1t),1);
a1t     = alpha1*ones(length(time2t),1);
a2t     = alpha0*ones(length(time3t),1);
plot([time1t;time2t;time3t],[a0t;a1t;a2t],'b-.','linewidth',line)

% marker minimum diluant : 
a3d     = alpha1*ones(length(time2d(25:27)),1);
a4d     = alpha1*ones(length(time2d(27:29)),1);
plot(time2d(27:29),a4d,'g--','linewidth',line)
plot(time2d(25:27),a3d,'k.','markersize',marker_size,'linewidth',line)
plot([time2d(29);time2d(29)],[1;1.07],'g--',time2d(29),1.07,'^','MarkerEdgeColor','g','MarkerFaceColor','g','markersize',marker_size)

% marker eco-optimal : 
plot([time2m(end);time2m(end)],[1;1.07],'r',time2m(end),1.07,'^','MarkerEdgeColor','r','MarkerFaceColor','r','markersize',marker_size)

% arrow for minimum time : 
% arrow([10;0.1],[4.179;0],1,'Ends','stop','FaceColor','k','EdgeColor','k','LineStyle','-','TipAngle',12,'Length',12,'linewidth',line)
arrow([8;0.2],[4.4;0],1,'Ends','stop','FaceColor','k','EdgeColor','k','LineStyle','-','TipAngle',12,'Length',12,'linewidth',line)
text(8.2, 0.22  , 'mintf', 'Color', 'k');

xlabel('time')
ylabel('alpha')
axis([0 18.9 -0.05 1.1])
legend('minimum diluant','economically optimal','minimum time',...
       'Location','SouthEast')
%fs=6; fz=11;
%exportfig(gcf,'RaEtCl_control.eps','width',fs,'fontmode','fixed',...
%         'fontsize',fz,'Color','rgb'); 
set(gcf, 'PaperUnits', 'centimeters') 
set(gcf, 'PaperPosition', [3 10 16 10])
print -depsc 'RaEtCl_control.eps'

function fun = model(t,x,u,param)
c1 = x(1); c2 = x(2); alpha = u;

q  = -29.19*log(c1) + 118.1;
R2 = 1 - (0.0073*c1 + 0.813);

dc1dt = param.K*c1^2*q*(1 - alpha);
dc2dt = c1*c2*param.K*q*(R2 - alpha);
daqdt = alpha*q;

fun = [dc1dt;dc2dt;daqdt];

function [value,isterminal,direction] = event1d(t,y,u,param)
% value      = param.b(4) - (-2*param.b(3)*y(1));
value      = y(1) - (param.clim - param.eps);
isterminal = 1;
direction  = 0;

function [value,isterminal,direction] = event1m(t,y,u,param)
value      = - param.wT*((param.b(4) + param.b(3)*y(1))*(param.b(1) +...
             param.b(2) + param.b(1)*log(y(1))) + param.b(3)*y(1)*...
             (param.b(2) + param.b(1)*log(y(1)))) - param.wD*...
             (param.b(2) + param.b(1)*log(y(1)))^2*(param.b(4) +...
             2*param.b(3)*y(1));
isterminal = 1;
direction  = 0;

function [value,isterminal,direction] = event1t(t,y,u,param)
% value      = param.b(2) - (-((param.b(1) + param.b(1)*log(y(1)))*(param.b(4) + param.b(3)*y(1)) + param.b(1)*param.b(3)*y(1)*log(y(1)))/(param.b(4) + 2*param.b(3)*y(1)));
value      = y(1) - param.copt;
isterminal = 1;
direction  = 0;

function [value,isterminal,direction] = event2d(t,y,u,param)
value      = (y(1)/y(2)) - (param.c1f/param.c2f);
isterminal = 1;
direction  = 0;

function [value,isterminal,direction] = event2m(t,y,u,param)
value      = (y(1)/y(2)) - (param.c1f/param.c2f);
isterminal = 1;
direction  = 0;

function [value,isterminal,direction] = event2t(t,y,u,param)
value      = y(2) - param.c2f;
isterminal = 1;
direction  = -1;

function [value,isterminal,direction] = event3d(t,y,u,param)
value      = y(1) - param.c1f;
isterminal = 1;
direction  = 0;

function [value,isterminal,direction] = event3m(t,y,u,param)
value      = y(1) - param.c1f;
isterminal = 1;
direction  = 0;

function [value,isterminal,direction] = event3t(t,y,u,param)
value      = y(1) - param.c1f;
isterminal = 1;
direction  = 0;