function states
%%% arrow.m needed
%%% http://www.mathworks.de/matlabcentral/fileexchange/278-arrow-m
  % initial and final conditions :
  x0 = [26; 16];           param.c10 = x0(1);  param.c20 = x0(2);   % [g/L]
  xf = [130;6.8];          param.c1f = xf(1);  param.c2f = xf(2);   % [g/L]
  param.V0 = 100;                                                   % [L]
  
  % parameters : 
  param.k    = 16.98;                                               % [m/h]                      
  param.A    = 1;                                                   % [m2]
  param.clim = 260;                                                 % [g/L]
  param.copt = param.clim/exp(1);                                   % [g/L]
  param.K    = 1/(param.c10*param.V0);
  param.P    = 0.95;
  param.eps  = 1e-4;
  param.N    = 10000;
  param.wT   = 3.9;
  param.wD   = 10;
  param.mulcopt = param.clim*exp((param.wT - (param.wT^2 + 4*param.k*param.A*param.wT*param.wD)^(1/2))/(2*param.A*param.k*param.wD));
  
  % control :
  alpha0   = 0;
  alpha1   = 1;
  alphainf = 1000;
  
  % Tolerances :
  AbsTol = 1e-13;
  RelTol = 1e-12;
  
  % integration step :
  step = 0.5;
  
  % figures options : 
  marker_size = 7;
  
  %% minimum time : 
  x01t      = [param.c10;param.c20;0];
  options1t = odeset('Events',@event1t,'AbsTol',AbsTol,'RelTol',RelTol);
  tspan1t   = [0:step:param.N];
  [time1t,state1t,te1t] = ode45(@model,tspan1t,x01t,options1t,alpha0,param);           
  
  options2t = odeset('Events',@event2,'AbsTol',AbsTol,'RelTol',RelTol);
  x02t      = state1t(end,:);
  tspan2t   = [time1t(end):step:param.N];
  [time2t,state2t,te2t] = ode45(@model,tspan2t,x02t,options2t,alpha1,param);
  
  %% minimum diluant :
  
  x01d      = [param.c10;param.c20;0];
  options1d = odeset('Events',@event1d,'AbsTol',AbsTol,'RelTol',RelTol);
  tspan1d   = [0:step:param.N];
  [time1d,state1d,te1d] = ode45(@model,tspan1d,x01d,options1d,alpha0,param);           
  
  options2d = odeset('Events',@event2,'AbsTol',AbsTol,'RelTol',RelTol);
  x02d      = state1d(end,:);
  tspan2d   = [time1d(end):step:param.N];
  [time2d,state2d,te2d] = ode45(@model,tspan2d,x02d,options2d,alpha1,param);
  
  %% multi-objective : 
  
  x01m      = [param.c10;param.c20;0];
  options1m = odeset('Events',@event1m,'AbsTol',AbsTol,'RelTol',RelTol);
  tspan1m   = [0:step:param.N];
  [time1m,state1m,te1m] = ode45(@model,tspan1m,x01m,options1m,alpha0,param);           
  
  options2m = odeset('Events',@event2,'AbsTol',AbsTol,'RelTol',RelTol);
  x02m      = state1m(end,:);
  tspan2m   = [time1m(end):step:param.N];
  [time2m,state2m,te2m] = ode45(@model,tspan2m,x02m,options2m,alpha1,param);
  
  %% almost minimum time :
  
  wTamt = 1;      % euro/h 
  wDamt = 0.1;    % euro/m3 
  param.coptamt = param.clim*exp((wTamt - (wTamt^2 + 4*param.k*param.A*wTamt*wDamt)^(1/2))/(2*param.A*param.k*wDamt));
  
  x01amt      = [param.c10;param.c20;0];
  options1amt = odeset('Events',@event1amt,'AbsTol',AbsTol,'RelTol',RelTol);
  tspan1amt   = [0:step:param.N];
  [time1amt,state1amt,te1amt] = ode45(@model,tspan1amt,x01amt,options1amt,alpha0,param);           
  
  options2amt = odeset('Events',@event2,'AbsTol',AbsTol,'RelTol',RelTol);
  x02amt      = state1amt(end,:);
  tspan2amt   = [time1amt(end):step:param.N];
  [time2amt,state2amt,te2amt] = ode45(@model,tspan2amt,x02amt,options2amt,alpha1,param);
  
  %% real minimum diluant :
  wTr = 0.01;      % euro/h 
  wDr = 1;         % euro/m3 
  param.coptrmd = param.clim*exp((wTr - (wTr^2 + 4*param.k*param.A*wTr*wDr)^(1/2))/(2*param.A*param.k*wDr));
  
  x01r      = [param.c10;param.c20;0];
  options1r = odeset('Events',@event1r,'AbsTol',AbsTol,'RelTol',RelTol);
  tspan1r   = [0:step:param.N];
  [time1r,state1r,te1r] = ode45(@model,tspan1r,x01r,options1r,alpha0,param);           
  
  options2r = odeset('Events',@event2,'AbsTol',AbsTol,'RelTol',RelTol);
  x02r      = state1r(end,:);
  tspan2r   = [time1r(end):step:param.N];
  [time2r,state2r,te2r] = ode45(@model,tspan2r,x02r,options2r,alpha1,param);
  
  %% figures :
  marker_sizes = 4;
  
  plot([state1d(:,1);state2d(:,1)],[state1d(:,2);state2d(:,2)],'r-','MarkerEdgeColor','r','MarkerFaceColor','r','markersize',marker_sizes)               % minimum diluant
  hold on
  plot([state1r(:,1);state2r(:,1)],[state1r(:,2);state2r(:,2)],'co-','MarkerEdgeColor','c','MarkerFaceColor','c','markersize',marker_sizes)              % realistic minimum diluant
  plot([state1m(:,1);state2m(:,1)],[state1m(:,2);state2m(:,2)],'kd-','MarkerEdgeColor','k','MarkerFaceColor','k','markersize',marker_sizes)              % economically optimal
  plot([state1amt(:,1);state2amt(:,1)],[state1amt(:,2);state2amt(:,2)],'gx-','MarkerEdgeColor','g','MarkerFaceColor','g','markersize',marker_sizes)      % almost minimum time
  plot([state1t(:,1);state2t(:,1)],[state1t(:,2);state2t(:,2)],'bs-','MarkerEdgeColor','b','MarkerFaceColor','b','markersize',marker_sizes)              % minimum time
  
  % purity : 
  c1 = [0:10:300];
  c2 = c1*((1 - param.P)/param.P);
  plot(c1,c2,'r--','MarkerEdgeColor','r','MarkerFaceColor','r')
  
  plot(param.c10,param.c20,'go','markersize',marker_size)
  arrtype = 13;
  arrow([state2t(end,1);state2t(end,2)],[state2t(end,1);state2t(end,2)],arrtype,'Ends','start','FaceColor','b','EdgeColor','b')
  arrow([state2m(end,1);state2m(end,2)],[state2m(end,1);state2m(end,2)],arrtype,'Ends','start','FaceColor','k','EdgeColor','k')
  arrow([state2d(end,1);state2d(end,2)],[state2d(end,1);state2d(end,2)],arrtype,'Ends','start','FaceColor','r','EdgeColor','r')
  arrow([state2r(end,1);state2r(end,2)],[state2r(end,1);state2r(end,2)],arrtype,'Ends','start','FaceColor','c','EdgeColor','c')
  arrow([state2amt(end,1);state2amt(end,2)],[state2amt(end,1);state2amt(end,2)],arrtype,'Ends','start','FaceColor','g','EdgeColor','g')
  xlabel('c1')
  ylabel('c2')
  legend('minimum diluant','realistic minimum diluant','economically optimal','almost minimum time','minimum time','P(c_1,c_2) = P_f','Location','SouthEast')
  axis([0 300 2 17])
  
  set(gcf, 'PaperUnits', 'centimeters')
  set(gcf, 'PaperPosition', [3 10 16 8])
  print -depsc 'soybean_state_purity.eps'

function fun = model(t,x,u,param)
  c1 = x(1); c2 = x(2); alpha = u;
  q = param.k*param.A*log(param.clim/c1);

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

  fun = [dc1dt;dc2dt;daqdt];
  
function [value,isterminal,direction] = event1t(t,y,u,param)
  value      = y(1) - param.copt;
  isterminal = 1; direction  = 1;  

function [value,isterminal,direction] = event2(t,y,u,param)
  value      = (y(1)/y(2)) - (param.P/(1 - param.P));
  isterminal = 1; direction  = 0;
  
function [value,isterminal,direction] = event1d(t,y,u,param)
  value      = y(1) - (param.clim*(1 - param.eps));
  isterminal = 1;  direction  = 1;
  
function [value,isterminal,direction] = event1m(t,y,u,param)
  value      = y(1) - param.mulcopt;
  isterminal = 1; direction  = 1;  
  
function [value,isterminal,direction] = event1amt(t,y,u,param)
  value      = y(1) - param.coptamt;
  isterminal = 1; direction  = 1;

function [value,isterminal,direction] = event1r(t,y,u,param)
  value      = y(1) - param.coptrmd;
  isterminal = 1; direction  = 1;  