function copt
%% 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    = 10000000000;
param.b    = [-29.19; 118.1; 0.0073; 0.813];
param.clim = 57.1630;
param.copt = 24.5001;
param.eps  = 1e-3;

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

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

% weight coefficients for calulating multi-objective 
% optimal time and % consumption
param.wT = 2;
wD       = 10;
%% Simulation :

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

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

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

% minimum diluant : 

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

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

% multi-objective

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

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

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

clear wD;
param.wT = 0.01;      % euro/h 
wDrd = 1;             % euro/m3 

% realistic minimum diluant : 

x01rd      = [param.c10;param.c20;0];
options1rd = odeset('Events',@event1m,'AbsTol',AbsTol,'RelTol',RelTol);
tspan1rd   = [0 param.N];
[time1rd,state1rd,te1rd] = ode45(@model,tspan1rd,x01rd,options1rd,alpha0,param,wDrd);           

options2rd = odeset('Events',@event2m,'AbsTol',AbsTol,'RelTol',RelTol);
x02rd      = state1rd(end,:);
tspan2rd   = [time1rd(end) param.N];
[time2rd,state2rd,te2rd] = ode45(@model,tspan2rd,x02rd,options2rd,alpha1,param,wDrd);

options3rd = odeset('Events',@event3m,'AbsTol',AbsTol,'RelTol',RelTol);
x03rd      = state2rd(end,:);
tspan3rd   = [time2rd(end) param.N];
[time3rd,state3rd,te3rd] = ode45(@model,tspan3rd,x03rd,options3rd,alphainf,param,wDrd);

% almost minimum time :
clear wD;
param.wT = 1;     % euro/h 
wDm1 = 0.1;       % euro/m3

x01amt      = [param.c10;param.c20;0];
options1amt = odeset('Events',@event1m,'AbsTol',AbsTol,'RelTol',RelTol);
tspan1amt   = [0 param.N];
[time1amt,state1amt,te1amt] = ode45(@model,tspan1amt,x01amt,options1amt,alpha0,param,wDm1);           

options2amt = odeset('Events',@event2t,'AbsTol',AbsTol,'RelTol',RelTol);
x02amt      = state1amt(end,:);
tspan2amt   = [time1amt(end) param.N];
[time2amt,state2amt,te2amt] = ode45(@model,tspan2amt,x02amt,options2amt,alpha1,param,wDm1);

options3amt = odeset('Events',@event3m,'AbsTol',AbsTol,'RelTol',RelTol);
x03amt      = state2amt(end,:);
tspan3amt   = [time2amt(end) param.N];
[time3amt,state3amt,te3amt] = ode45(@model,tspan3amt,x03amt,options3amt,alpha0,param,wDm1);

%% Pareto-Front
clear wD;
% weigth coeficients : 
param.wT = 1;
wD       = logspace(-6,6,100);

% memory allocation :
tff      = zeros(length(wD),1);
Vw       = zeros(length(wD),1);
c1opt    = zeros(length(wD),1);

% multi-objective :
for i = 1 : length(wD)

   x01mm      = [param.c10;param.c20;0];
   options1mm = odeset('Events',@event1m,'AbsTol',AbsTol,'RelTol',RelTol);
   tspan1mm   = [0 param.N];
   [time1mm,state1mm,te1mm] = ode45(@model,tspan1mm,x01mm,options1mm,alpha0,param,wD(i));           

   if (state1mm(end,1) < param.c1f)
   
     options2mm = odeset('Events',@event21m,'AbsTol',AbsTol,'RelTol',RelTol);
     x02mm      = state1mm(end,:);
     tspan2mm   = [time1mm(end) param.N];
     [time2mm,state2mm,te2mm] = ode45(@model,tspan2mm,x02mm,options2mm,alpha1,param,wD(i));

     options3mm = odeset('Events',@event3m,'AbsTol',AbsTol,'RelTol',RelTol);
     x03mm      = state2mm(end,:);
     tspan3mm   = [time2mm(end) param.N];
     [time3mm,state3mm,te3mm] = ode45(@model,tspan3mm,x03mm,options3mm,alpha0,param,wD(i));
   
   elseif (state1mm(end,1) > param.c1f)
       
     options2mm = odeset('Events',@event22m,'AbsTol',AbsTol,'RelTol',RelTol);
     x02mm      = state1mm(end,:);
     tspan2mm   = [time1mm(end) param.N];
     [time2mm,state2mm,te2mm] = ode45(@model,tspan2mm,x02mm,options2mm,alpha1,param,wD(i));
     
     options3mm = odeset('Events',@event3m,'AbsTol',AbsTol,'RelTol',RelTol);
     x03mm      = state2mm(end,:);
     tspan3mm   = [time2mm(end) param.N];
     [time3mm,state3mm,te3mm] = ode45(@model,tspan3mm,x03mm,options3mm,alphainf,param,wD(i));
     
   end
     
   tff(i)    = te3mm;
   Vw(i)     = state3mm(end,3);
   c1opt(i)  = state1mm(end,1);
end
%% Figures :

line        = 1.5;
marker_size = 9;

figure
plot(te3t, state3t(end,3),'bs','markersize',marker_size)                   % minimum time
hold on 
plot(te3rd, state3rd(end,3),'bo','markersize',marker_size)                 % realistic minimum diluant
plot(te3m,state3m(end,3),'bd','markersize',marker_size)                    % multi-objective
plot(te3amt,state3amt(end,3),'bh','markersize',marker_size)                % almost minimum time
plot(te3t,state3d(end,3),'b*','linewidth',line)                            % utopia point
plot([te3t;te3t],[Vw(1);Vw(end)],'b--','linewidth',line)                   % J = tf
plot([te3t;te3d],[state3d(end,3);state3d(end,3)],'b--','linewidth',line)   % J = Vw
plot(tff, Vw,'linewidth',line)
xlabel('tf')
ylabel('Vw')
legend('minimum time','realistic minimum diluant','economically optimal','almost minimum time','utopia point')
axis([0 60 22 80])
% fs=6; fz=11;
% exportfig(gcf,'RaEtGl_PF.eps','width',fs,'fontmode','fixed','fontsize',fz,'Color','rgb');
set(gcf, 'PaperUnits', 'centimeters') 
set(gcf, 'PaperPosition', [3 10 11 8.5])
print -depsc 'RaEtGl_PF.eps'

figure
cg1 = param.clim*ones(size(wD));
cg2 = param.copt*ones(size(wD));
semilogx(param.wT./wD,cg1,'r--',param.wT./wD,c1opt,'b',param.wT./wD,cg2,'g-.','linewidth',line);
%hold on
%semilogx([param.wT./wD(55);param.wT./wD(32)],[57.16;24.5],'k','linewidth',line)
%hold on
%semilogx([wD(1);param.wT./wD(55)],[57.16;57.16],'k','linewidth',line)
%semilogx([param.wT./wD(32);wD(end)],[24.5;24.5],'k','linewidth',line)
xlabel('wTwD')
ylabel('c1')
legend('S_D(c_1) = 0','c_1^*','S_T(c_1) = 0','Location','West')
%fs=6; fz=11;
%exportfig(gcf,'RaEtGl_c1opt_final.eps','width',fs,'fontmode','fixed','fontsize',fz,'Color','rgb');
set(gcf, 'PaperUnits', 'centimeters') 
set(gcf, 'PaperPosition', [3 10 16 10])
print -depsc 'RaEtGl_c1opt.eps'


% %% Production cost for 1000 kg :
% wT = 2;
% wD = 10;
% 
% Jt = (wT*time3t(end)   + wD*state3t(end,3))/0.48*10
% Jm = (wT*time3m(end)   + wD*state3m(end,3))/0.48*10
% Jd = Inf;
% Jr = (wT*time3rd(end)  + wD*state3rd(end,3))/0.48*10
% Ja = (wT*time3amt(end) + wD*state3amt(end,3))/0.48*10
% Jccvd = (wT*4.22 + wD*93.31)/0.48*10
% % Jvvd  = (wT*4.3353 + wD*90.9793)/0.48*10






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

q  = param.b(1)*log(c1) + param.b(2);
R2 = 1 - (param.b(3)*c1 + param.b(4));

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,wD)
% 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,wD)
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)))) - 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,wD)
% 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,wD)
value      = (y(1)/y(2)) - (param.c1f/param.c2f);
isterminal = 1;
direction  = 0;

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

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

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

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

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

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

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