function table
% c1 = macro-solute (proteins)
% c2 = micro-solute (inpurities)

%% Type of problem: (Parameter "N")
% 1 --> minimum time
% 2 --> minimum diluant
% 3 --> multi-objective optimal

%% Initial/final conditions + paramters : 
x0 = [26; 16];     c10 = x0(1); c20 = x0(2);  % [g/L]
V0 = 100;                                     % [L]                  

% parameters : 
k    = 16.98;                                 % [m/h]                      
A    = 1;                                     % [m2]
clim = 260;                                   % [g/L]
P    = 0.95;                                  % purity [%]

% weight coefficients : 
wT      = 1;                                  % -
wD      = logspace(-6,6,100);                 % -

% memort allocation :
copt = zeros(length(wD),1);
tf   = zeros(length(wD),1);
Vw   = zeros(length(wD),1);

%% Calculating optimal time and water consumption for different values of weight coefficients
nn = length(wD) ;
for i=1:nn
  copt(i) = clim*exp((wT - (wT^2 + 4*k*A*wT*wD(i))^(1/2))/(2*A*k*wD(i)));
  [tf(i),Vw(i),c11,c12,c22,dt1,dt2]=find_tf_vw(wT,wD(i),c10,c20,V0,clim,k,A,P,3);
  str1 = fprintf('Step: %d/%d, wTwd: %d \n', i,nn,wT/wD(i));
end

%% Calulation of minimum time, dilunat, ecocnomically optimal, almost min time and real min diluant :

% calculating time and water consumption (minimum time problem): 
wTt = 1;
wDt = 0;
[tft,Vwt,c11t,c12t,c22t,dt1t,dt2t]=find_tf_vw(wTt,wDt,c10,c20,V0,clim,k,A,P,1);

% calculating time and water consumption (minimum diluant problem): 
wTd = 0;
wDd = 1;
[tfd,Vwd,c11d,c12d,c22d,dt1d,dt2d]=find_tf_vw(wTd,wDd,c10,c20,V0,clim,k,A,P,2);

% Calculating multi-objective optimal point :
wTm = 3.9;     % euro/h 
wDm = 10;      % euro/m3 
[tfm,Vwm,c11m,c12m,c22m,dt1m,dt2m]=find_tf_vw(wTm,wDm,c10,c20,V0,clim,k,A,P,3);

% Calculating almost minimum time :
wTm1 = 1;         % euro/h 
wDm1 = 0.1;       % euro/m3 
[tfm1,Vwm1,c11m1,c12m1,c22m1,dt1m1,dt2m1]=find_tf_vw(wTm1,wDm1,c10,c20,V0,clim,k,A,P,3);

% Calculating real minimum diluant :
wTd1 = 0.01;      % euro/h 
wDd1 = 1;         % euro/m3 
[tfd1,Vwd1,c11d1,c12d1,c22d1,dt1d1,dt2d1]=find_tf_vw(wTd1,wDd1,c10,c20,V0,clim,k,A,P,3);

c1=logspace(log10(10),log10(300),60);
wT=1; wD=eps;
Sc2t=clim*exp((wT - (wT^2 + 4*k*A*wT*wD)^(1/2))/(2*A*k*wD));

wT=0; wD=1;
Sc2d=clim*exp((wT - (wT^2 + 4*k*A*wT*wD)^(1/2))/(2*A*k*wD));

wT=1; wD=1;
Sc2m=clim*exp((wT - (wT^2 + 4*k*A*wT*wD)^(1/2))/(2*A*k*wD));

%% Pareto - front

marker_size = 7;
figure
plot(tft, Vwt,'bs','markersize',marker_size)                               % minimum time
hold on
plot(tfd1,Vwd1,'bo','markersize',marker_size)                              % minimum diluant
plot(tfm,Vwm,'bd','markersize',marker_size)                                % multi-objective optimal
plot(tfm1, Vwm1,'bx','markersize',marker_size)                             % almost min time 
plot(tft,Vwd,'b*','markersize',marker_size)                                % utopia point
plot(tf,Vw)
plot([tft;tft],[Vw(1);Vw(end)],'b--')                                      % J = tf
plot([tf(1);tf(end)],[Vwd;Vwd],'b--')                                      % J = Vw
legend('minimum time','realistic minimum diluant','economically optimal','almost minimum time','utopia point')
axis([4 12 0 32.4])
xlabel('tf')
ylabel('Vw')
%fs=6; fz=11;
%exportfig(gcf,'soybean_free.eps','width',fs,'fontmode','fixed','fontsize',fz,'Color','rgb');
set(gcf, 'PaperUnits', 'centimeters')
set(gcf, 'PaperPosition', [3 10 11 8.25])
print -depsc 'soybean_free.eps'

%% Results : 
clc
str2 = fprintf('Minimum time problem:    tf = %d [h], Vw = %d [L], J = %d [euro]\n',tft,Vwt,(wTm*tft + wDm*Vwt)/2.6*1000);
str3 = fprintf('Minimum diluant problem: tf = %d [h], Vw = %d [L], J = %d [euro]\n',tfd,Vwd,(wTm*tfd + wDm*Vwd)/2.6*1000);
str4 = fprintf('Multi-objective optimal: tf = %d [h], Vw = %d [L], J = %d [euro]\n',tfm,Vwm,(wTm*tfm + wDm*Vwm)/2.6*1000);
str5 = fprintf('Almost minimum time :    tf = %d [h], Vw = %d [L], J = %d [euro]\n',tfm1,Vwm1,(wTm*tfm1 + wDm*Vwm1)/2.6*1000);
str6 = fprintf('Real minimum diluant :   tf = %d [h], Vw = %d [L], J = %d [euro]\n',tfd1,Vwd1,(wTm*tfd1 + wDm*Vwd1)/2.6*1000);

function [tf,Vw,c11,c12,c22,delt1,delt2]=find_tf_vw(wT,wD,c10,c20,V0,clim,k,A,P,N)

if (N == 1)
  copt = clim/exp(1);
elseif (N == 2)
  copt = clim;
else
  copt = clim*exp((wT - (wT^2 + 4*k*A*wT*wD)^(1/2))/(2*A*k*wD));
end

% ---------------------------------------------------------------------------------

     % 1st phase (c2 = const, c10 -> c11) pure ultrafiltration (alpha = 0)
     % until c1 = copt
     
     c11   = copt;
     delt1 = ((c10*V0)/(k*A*clim))*(mfun('Li',clim/c10) - mfun('Li',clim/c11));
     if (N==2)
       delt1 = Inf;
     end
     Vw1   = 0;

% ---------------------------------------------------------------------------------

     % 2nd phase (c1 = const, c12 -> c11) diafiltration
     c12 = copt;

     % CVD step until : P = c1/(c1 + c2)
     c22   = c12*(1 - P)/P;
       if (N ~=2) 
	 delt2 = ((c10*V0)/(k*A*c12))*((log((c20*P)/(c12*(1 - P))))/(log(clim/c12))); 
	 Vw2   = (k*A*log(clim/c12))*delt2;
       else
	 delt2 = Inf;
	 Vw2 = (c10*V0/clim)*(log((c20*P)/(clim*(1 - P))));
       end

% ---------------------------------------------------------------------------------

     % final results : 
     tf = delt1 + delt2;
     Vw = Vw1 + Vw2;
%     [wT, wD, delt1 delt2 Vw1 Vw2]
