From 3d20a55f41046121d6ae728d058676af2ad1acf8 Mon Sep 17 00:00:00 2001
From: Oliver <oliver-mx@gmx.de>
Date: Fri, 9 Aug 2024 15:09:41 +0200
Subject: [PATCH] Parallel Matlab code

---
 DATA_case1.mat  | Bin 0 -> 2768 bytes
 DATA_case2.mat  | Bin 834 -> 826 bytes
 ODEsystem.m     |   4 +-
 Pareto_1_data.m | 125 +++++++++++++
 Pareto_2_data.m | 125 +++++++++++++
 Sim_0a_data.m   |   8 +-
 Sim_0b_data.m   |   8 +-
 Skript.asv      | 137 +++++++++++---
 Skript.m        | 144 +++++++++++++--
 fun_1.asv       | 455 ++++++++++------------------------------------
 fun_1.m         | 471 +++++++++---------------------------------------
 nonlcon.m       |   9 +-
 12 files changed, 688 insertions(+), 798 deletions(-)
 create mode 100644 DATA_case1.mat
 create mode 100644 Pareto_1_data.m
 create mode 100644 Pareto_2_data.m

diff --git a/DATA_case1.mat b/DATA_case1.mat
new file mode 100644
index 0000000000000000000000000000000000000000..004ec10e0d4293447c5ef96ff1b845873931e076
GIT binary patch
literal 2768
zcmeZu4DoSvQZUssQ1EpO(M`+DN!3vZ$Vn_o%P-2cQV4Jk_w+L}(NS<NN=+<DO;O0t
zvr=#?%2aSHO;=E`R4_EPGBB_*G*vJ%Ffvgf5_EGiKmensFFykV!wMz_hKf15XRAGe
z3q{yJSbJ{|IvtuNm3lMnl8Wv}!N|JHm$JU!yT3OjFwAL^Y9_O;arJ`}PF_NPoqajC
zJ~TOc;n2Jb-A7GIX7)$T+&t&m7NaJiz?j$9A0K!)=e*^gx^%h2250?v9y_h>`5FGd
z<#YXIQDKi@r`@*Ml4*Qzx0cL!zdCyg|Hc`e<ws;I&mP*Vr(jdsbu>9P+{0bJ`{%_r
zdHJa_dD6#^Dcfx>Fsj=r6mfc?j+2qOQlm59>X;v^<{!Sk<#Ff9nKFB}tL1k^o;_T*
zcU^q&@q;;cKNh-b^;c`Hzxb<+?=5R(Uo@{xQSYndOY^_X-IhHuaZ-Ct;n5#=8RV|^
zaxVKM(X^^A`^|?uIliyWE8o>LEqO30`}|g$Fwyue&u4Kwlln6E`wRQzb?Rpv7qd)T
z$`$Fa+t1cpSHGseLGA6XS&!o6vR>{!#l@jv*1mT|nPS1m?5XQch$ZMp2K`9ZKX`V|
zb>WAuOg3y?`<5>e_9$KY?6f~ihRycc`xksBg#Y+?`KuzMqi_3>Wp+(mkHtGzhq&Ck
zdUTsqC&SeR$7_QZaJc#}OD$N~(765Go0%1YyB3$3Khl2?vO>{PDfY;Z-w(dDPh=8n
z@oV36v17(F$@oPA3AIyV>y^A)(xbvlk0~+k{=C;N^=Vgx-VBrS6ec6KOHpkrXBy1^
z6vBP*=c?M<%U3Sm*{eSL3rpa!6PqS(R;s(|S@ug_C1BpS+eRMWJ9j&|6@{|gV)6f9
z`0A{%(Rx4W{v#6Avpa9;NBmUWIjP;iHaTbMy#N1h2AH2uPCe7hWL%oR-QWMTO8eo$
z{<*qWir#O_-|c0oGSy~kf6#9JB`x3AKDv-JckSlD-8WB2r|NbunZo?&e)TWzyQX?N
zRx+Wl)mP{rooM<@!F1Z*IigP&d<&a<I&SrizWXOvDAt40s}=(&y*^G(NXTGFGE;cg
za7p0+*HdH0jcfwOMjH%VB#yB+{9xp{!p`tzBU2_whaD?Ohwl|`gDcr03?KHVpWdp!
zJ}+$6#Z@OBw@tA*m$W-@@~ySoyvv1`2HjEl^5D|jcPjl}s?jgaPH9wzEQ;3H7rWLy
z?sS#TCY?pIG`Lni$=db!_3Za&ex5s%f9}jo%g=!vO&nKMI%axL=y?D0_QfCj|4o~I
zriE*v*AE?Kj<qw*v)Qhi<kcU0_*!B1%9jlVk24&~LzZuL;ry^CY`w>>Ta&-;J$*d$
zxwHIPo5Q6SV$^0=oBWiyvp>S-P23E=Ss#BjzLAevv_EWZ*&;rvQe*oW!V}qw78V@6
za3?1qGJo%;1nsh~du3-uuY0{<uGrSQm%ryRy!U>-AUx8p=0x>#fhqlWGq?7JHM?59
znDb(FlFehWno7ga?Ius+ExujewZ(p3qPU=2@Ab(q7gVe(bqw{=d$UzLt3D}$X*uto
z9QR`zs;#%~zgp<QezfhwG=ur~n1j>%G#I8j>Rp_B^3vy7rN6gc{lTYVBviHLt8amJ
zgYou~e7$JJwl^~8R;gIOi3@+?dvRvUmv3VCk3@fB2*~<hcl77-1+vG&IkqQ0-&68t
zg?x9<jg^bzYj5@+m}&RN=hG&pBd%gnJ^MFs+eIg)9adMK8Ty1VJJ2v?>%<v-epek<
zo}RpUU9x$~tKXLYBqUQ;znB^jJO6;{_uRL$Cucc0XK!!NIiyz^I;VQUgdU;P9)GsX
zREDg{^X^=-@>njuZ)I)MU4!R-2R|?ybw!<N5f%}gW3<mFW#W<rTWW;;UmZNR;_j+B
zx18Mk^!QC{ER!ajcAqPL)AZ1NJ!#%$SFfMjRiml7{4=wpDoe~2tLmLHdqU+l@BTY`
zfrBkW-tNTbsz$|~SC=pFd2hJ1{&a9g5pV2PM}_A2b26NrACKv^i-*q4^mqIgd)q>v
zM_VB_$Nb^c%RF*_PCT$r-f((rrm6XDC&$vJ2)9$yybQ$m@=m(G{lKzy*K0DR%8JjQ
zzrRq{)A6=GOL0u_Duo$4A5JJ(BFc4m-<zJA4%x%eTi=)~vZ}HNc`jO(rO|bcm;b`#
z!+k0L&olGhczfz=P=r|Eg&MZGCl_CLi1?INwxHIUM{sT?@7^1|<!^3^-}ZA<ck0?K
znYZ=9guvQ)dzUfWam2mlRJr!hcJHHQ(=B`*FP#uIY~w7w5cKN#mhTSg7RTSs`|$Sg
z-%BP($}758aIo5mMTcB<{*=_Zs{T^(4b$_<OH*~rCjDA9>FMj0Cz5`>Sn?p4nI}I@
z-}C=+?p5+DtJ%bo%_SF}iIZ@3HS1B__R3WA&Y8lfRbguyKCQe`5XDu-zTVY9vA-)}
z{j=3|)7@?ypQ~<v^$KtM-POrwz9nAYQt;b*s-sqw{>{o+H!pKA_Z6Pn*dceKru@E1
z@=yL^k>b~nGEA<Nn`uv%zRF^*XX+;Z^A7X&OeXeDv$9#jV*5|$WUfd*l9qD7gRM|4
z^-7Ph*8|rt-=3WQFmGS<0j=}j7b@Cco#<D#k*7rac8<;SFHN7kl4l2MmF<q{o!Y0K
zsnIADSbTlCq>86BPtW7(hjUiCO%P<?T`=!P+MbR!K8Y3n+5TqDTjz?+w_VQ9XSzmI
zBV=Ro?t=2{7c9JMzy13eacrrWfN<BwXOC+itl`*w=+vSA7p(4I;#4lTN}76OTR`}>
zrb$!ELO<=}UsU}(Gp6*QUvWdVjn&+~p1oeyMz^C917C#J%h%cW7#6<_Rn1PjtC2q;
zCVf%&vDGabJ`_mLNeQ3z+~y7Y?Ui3OXY*`K?pb<8r0k_q_-c*n#)UFFf^Kh`{{FY)
ztWt~5^PaBXoOiIsdxFo)#YY`$uD@e1y0HD!j<dTi@M~Uc@jO4__?qARQ>Lvs0xCtd
z;pOfR#v?fnXBd|{Hi({LzbNH!=EQ+R0+|x>2Ut?(CrH>aFub>Bjs<D)hiiEOE~*)l
z9N9P$SdzjN^iH1fHD&o6DN!@`Z?ojWS3zFgLf%H7KmTaR;aB;)C*bdz06Xa`I&Ld0
zI`kP6@1;6TYEER})Z{n`GByy0v1Jfrd-xeXkz(*h#-Mkc)gXgy89@~VsA};@NKHs!
zz9i#t<|G5-ngBax0qG4+83GASpCl!0YaMIs6y_ahx^TY5)t<-sNaTT{o&eK*=l&fu
zD7Kk!MzCd7ddp><gslp8`U1OCpDk$<R}~XkB6jq+^Ved9m+m*%A8uZ}lX*#8qu}=a
z_6t^~H+3<tO-Nu%Vh`Bi_>oItAHxy-1rkg2%#GF6H_uR9x<zB7im`paqVDz<M)8D0
zn#~M75kGn&j?BDp#yR|Whbxc#1LsA$5&{p8vWXrSJwGSlV;l2wsrTnL6{&>VJ0CtX
zbJE4vG7i7O4IcExtFfK=;ku1KdRIz<gYWb$SI#>h*lSS8mh!jf%h&guOusIEHq`NN
z{Ni-k+Ggft+4~lM4mY{lZ`#kY;7<?F>N_#3A0<ur$H(^PX%x$(IyDXXoiU#k3~dUg
NW-Bl->^fq28vp}_4R`<m

literal 0
HcmV?d00001

diff --git a/DATA_case2.mat b/DATA_case2.mat
index f14208f56b7a771b6bfd8753081021c22ffb5451..e66e8aa3eabf7b54bb1562ff923c4f07abfb300f 100644
GIT binary patch
delta 156
zcmX@awu^0oiG*8Grh;Q>x`KkGf}xR>k(rf&p@NZtk;%kB<%tQbOn(_CmX)S!GcYhz
z%z2!gknn@?NRGo9#-)x8qNms|N;#Z4ao~tR<`D)JNeMp#9(e|acon9#{0s~Xc(oi7
i$P|%3z>+dQLBfuK;R7FY97v1*<c*Btn+q8)G6Dck*fNO#

delta 155
zcmdnRc8G0)i9|?7se)r^x`Kj*f}w$xv4NF=k%Ezdk;%kB<%tQbOsq^3%Sy9!7#J8T
z<~&YLNZ7$}B*)<l<5I^4(NpP+NxTYaDG5m(CP@sFczC!a*z6e?c1&f8;b&lAFoNsh
zNJvXaV7?^7BOt{Rv9Lx=!H)4rk3-m!z#0bO32P6WG2n7LsBwUqVc8<4&6|B0-!cLK
DqR1}(

diff --git a/ODEsystem.m b/ODEsystem.m
index 77862a4..b4c7306 100644
--- a/ODEsystem.m
+++ b/ODEsystem.m
@@ -31,7 +31,7 @@ swro_del_c = J_p(1)./(J_p(1) + 1) - J_p(3)./(J_p(3) + J_p(4));     % delta c_sal
 J_cross = ((J_p(5) - J_p(6)) - sigma.*(swro_p_osm_d - swro_p_osm_f)) ./   (1 + p_r*swro_alpha*swro_KK*sigma.*(swro_p_osm_d - swro_p_osm_f));
 if version(4)==0; J_cross =((J_p(5) - J_p(6)) - sigma.*(swro_p_osm_d - swro_p_osm_f)); end
 J_sin = swro_beta*swro_del_c;
-if version(7)==1 %% version with ICP and ECP
+if version(7)==1 && version(4)>0
     J_cross = ((J_p(5)-J_p(6))-sigma.*(swro_p_osm_d - swro_p_osm_f) + swro_beta*(J_p(5)-J_p(6))*(swro_KK +1/swro_KD + 1/swro_KF)) ./   (1 + swro_alpha*swro_beta*(1/swro_KD + 1/swro_KF + swro_KK) + p_r*swro_alpha*sigma.*(-swro_p_osm_d/swro_KD -swro_p_osm_f*(swro_KK + 1/swro_KF) ));
     J_sin = swro_beta*( swro_del_c-J_cross*( J_p(1)./(J_p(1) + 1)./swro_KD +  J_p(3)./(J_p(3) + J_p(4))*(swro_KK+1/swro_KF)))/(1+swro_alpha*swro_beta*(swro_KK+1/swro_KF+1/swro_KD));
 end
@@ -66,7 +66,7 @@ del_c = J_p(1+6)./(J_p(1+6) + 1) - J_p(3+6)./(J_p(3+6) + J_p(4+6)); % delta c_sa
 Q_cross = ((p_osm_d - p_osm_f) - (J_p(5+6) - J_p(6+6)) - (J_p(5+6) - J_p(6+6))*beta*KK*alpha*p_r)./(1 + KK*alpha*p_r*(beta +p_osm_f));
 if version(5)==0; Q_cross = ((p_osm_d - p_osm_f) - (J_p(5+6) - J_p(6+6))) ; end 
 Q_sin = beta*del_c;
-if version(8)==1 %% version with ICP and ECP
+if version(8)==1 && version(4)>0
     Q_cross=((p_osm_d - p_osm_f) - (J_p(5+6) - J_p(6+6)) - (J_p(5+6) - J_p(6+6))*beta*(KK+1/KF+1/KD)*p_r*alpha)./(1 +alpha* beta*(KK+1/KF+1/KD) + alpha*p_r*(p_osm_f*(KK+1/KF) + p_osm_d/KD));
     Q_sin=beta*( del_c-Q_cross*( J_p(1+6)./(J_p(1+6) + 1)./KD +  J_p(3+6)./(J_p(3+6) + J_p(4+6))*(KK+1/KF)))/(1+alpha*beta*(KK+1/KF+1/KD));
 end
diff --git a/Pareto_1_data.m b/Pareto_1_data.m
new file mode 100644
index 0000000..1c37970
--- /dev/null
+++ b/Pareto_1_data.m
@@ -0,0 +1,125 @@
+function [H,Z,swro_Z,ro_water,ro_salt,Mw,Ms,Rw,T0,eta,sigma,p_r,rho_r,C_r,swro_L,swro_alpha,swro_R,swro_KK,swro_x_r,swro_b1,swro_b2,J_r,swro_gamma,swro_gamma2,swro_W_r,L,alpha,R,KK,x_r,b1,b2,Q_r,gamma,gamma2,W_r,cE,pE,rho_E,J_sf_0,J_wf_0,Pd_0,Pd_L,Pf_L,Q_sf_0,pd_0,pf_0,pd_L,pf_L,HP_eff,LP_eff,T_eff,V_m,ERD_eff,ERD_fric,A_ERD,eta_ERD,mix_density,pw,pe,swro_beta_fix,beta_fix,mixer_ERD,version,fig,swro_KF,swro_KD,KF,KD]= Pareto_1_data(input1, input2)
+%%  Pareto_2_data(input)
+%
+%   Data for Berechung der roten Pareto front
+%
+%   only SWRO
+%
+%   option_data = 1
+%
+%   Input:
+%       input1        -   x
+%       input2        -   x
+
+%% model versions
+version=zeros(1,10);
+% version(1)=0 if co-current, 1 otherwise
+version(2)=1;  % 0 = SWRO beta fixed
+version(3)=0;  % 0 = PRO beta fixed
+version(4)=1;  % 0 = ideal SWRO
+version(5)=1;  % 0 = ideal PRO
+% ERD
+version(6)=0;
+% version(6) = 0 --> only SWRO (no ERD)
+% version(6) = 1 --> only SWRO (with ERD)
+% version(6) = 2 --> only PRO
+% version(6) = 3 --> SWRO-PRO hybrid system (with one ERD)
+% version(6) = 4 --> SWRO-PRO hybrid system (with two ERDs) 
+version(7)=1; % 1 = ICP and ECP for SWRO
+version(8)=1; % 1 = ICP and ECP for PRO
+
+%% Membrane unit properties
+H = 1e-3;           % height of the membrane [m]
+ro_water = 1000;    % mass density of water  [kg/m^3] 
+ro_salt = 2165;     % mass density of salt [kg/m^3] 
+Mw = 18;            % molecular weight of water [kg/kmol]
+Ms = 58.44;         % molecular weight of salt  [kg/kmol]
+Rw = 462;           % gas constant of water  [J/(kg K)] 
+T0 = 297;           % temperature [K] 
+eta = 1.3e-3;       % seawater viscosity [kg/(m s)]
+p_r = 1e5;          % pressure [Pa]=[kg/ms^2]
+rho_r = 1e3;        % density [kg/m^3]
+C_r = 1;            % salt concentration [%]
+
+%% SWRO
+swro_Z=1;               % width of SWRO membrane [m]
+swro_L=4;               % length of SWRO membrane [m]
+if version(6)== 2; swro_L=1; end
+swro_alpha = 5.0815e-9; % SWRO water permeablity co-efficient [s/m]
+swro_R =0.96;           % SWRO salt rejection rate 
+swro_KK = 1e-2;         % SWRO ICP mass transfer coefficient
+swro_KD = 1/swro_KK;    % SWRO ECP draw side mass transfer coefficient
+swro_KF = 1/swro_KK;    % SWRO ECP fresh side mass transfer coefficient
+swro_x_r=swro_L^2;      % x_r=swro_L^2 since x = linspace(0,1,n) (if x = linspace(0,swro_L,n) then x_r=swro_L;)     
+swro_b1 = H/swro_x_r;                           % H/swro_L ratio
+swro_b2 = swro_Z/swro_x_r;                      % Z/swro_L ratio
+J_r = sqrt(H^3/swro_x_r*p_r*rho_r);             % flux [kg/s^2]
+swro_gamma= swro_x_r * p_r *  swro_alpha /J_r;  % SWRO scaling factor - mass balance
+swro_gamma2= J_r^2./(swro_x_r^2 * p_r * rho_r); % SWRO scaling factor - momentum balance
+swro_W_r=J_r*p_r/rho_r;                         % net work [W/m^2]
+sigma=0.999 ;                                   % Rejection coefficient
+swro_beta_fix=4.43e-4/J_r*swro_x_r;             % value for fixed SWRO beta [kg/sm^2]
+ 
+%% PRO
+Z=1;                    % width of the PRO membrane [m]
+L=1.5;                  % length of the PRO membrane [m]      
+alpha = 5.47e-9;        % water permeablity co-efficient [s/m]
+R =0.94;                % salt rejection rate [1]
+KK=7.13e2;              % mass transfer coefficient [sm^2/kg]
+KD=1/KK;                % PRO ECP draw side mass transfer coefficient [sm^2/kg]
+KF=1/KK;                % PRO ECP fresh side mass transfer coefficient [sm^2/kg]
+x_r=L^2;                % x_r=L^2 since x = linspace(0,1,n) (if x = linspace(0,L,n) then x_r=L;)              
+b1 = H/x_r;             % H/L ratio
+b2 = Z/x_r;                           % Z/L ratio
+Q_r = sqrt(H^3/x_r*p_r*rho_r);        % flux [kg/s^2]
+gamma= x_r * p_r * alpha /Q_r ;       % PRO scaling factor - mass balance
+gamma2= Q_r^2./(x_r^2 * p_r * rho_r); % PRO scaling factor - momentum balance
+W_r=Q_r*p_r/rho_r;                    % net work [W/m^2]
+beta_fix =1.71e-4/Q_r*x_r;            % value for fixed PRO beta [kg/sm^2]
+
+%% Sea Water
+cE= 35/983/C_r;                            % salt concentration in seawater
+pE= 1e5/p_r;                               % external pressure
+rho_E=(cE + 1)./(ro_water*cE/ro_salt + 1); % density of incomming seawater
+
+%% SWRO operating conditions
+J_sf_0 = 0;              % salt flux in fresh side at 0 
+J_wf_0 = 0;              % water flux in fresh side at L               
+Pf_L = pE;               % pressure fresh side at L
+
+Pd_0 =  input1(1);       % pressure draw side at 0
+Pd_L =  input1(2);       % pressure draw side at L (not need in the hybrid system)
+
+%% PRO operation conditions 
+Q_sf_0 = 0;              % salt flux in fresh side at 0
+pf_L = pE;               % pressure of fresh side at L
+
+pd_0 = 12e5/p_r;         % pressure draw side at 0
+pd_L = 14.6e5/p_r;       % pressure of fresh side at 0
+pf_0 = 1.3e5/p_r;        % pressure draw side at L 
+
+%% ERD/Turbine/Pump parameters
+T_eff  = .95;               % turbine efficiency
+HP_eff = .9;            	% high pressure pump efficiency
+LP_eff = .95;               % low pressure pump efficiency
+V_m=0.052;                  % Volumetric mixing
+ERD_eff=.96;                % ERD unit pressure efficiency
+ERD_fric=5e-04;             % ERD friction coefficient
+A_ERD=H*swro_Z;             % cross sectional area of ERD inflows/outflows
+eta_ERD=0.01;               % leak of high pressure brine
+mix_density=997/rho_r;  	% density of mixture in ERD
+pw=2.4;                     % water price [$/m^3]
+pe=0.14;                    % electricity price [$/kWh]
+mixer_ERD=1;                % PRO Draw outlet mixer adjustment (only if 2nd ERDs) (mixer_ERD=1 --> all flow to ERD2 no turbine needed)    
+
+%% display figures
+%fig=[1,1,1,0]; % f(i)=1 --> figure i will be displayed
+fig=[0,0,0,0]; % f(i)=1 --> figure i will be displayed
+
+%% model specific changes:
+% co-current
+if pd_0 > pd_L; version(1)=0; else; version(1)=1; end            
+% no SWRO needed (--> use trivial data)
+if version(6) == 2; version(4)=0; Pd_0=40; Pd_L=35; end
+
+end
\ No newline at end of file
diff --git a/Pareto_2_data.m b/Pareto_2_data.m
new file mode 100644
index 0000000..2f910e7
--- /dev/null
+++ b/Pareto_2_data.m
@@ -0,0 +1,125 @@
+function [H,Z,swro_Z,ro_water,ro_salt,Mw,Ms,Rw,T0,eta,sigma,p_r,rho_r,C_r,swro_L,swro_alpha,swro_R,swro_KK,swro_x_r,swro_b1,swro_b2,J_r,swro_gamma,swro_gamma2,swro_W_r,L,alpha,R,KK,x_r,b1,b2,Q_r,gamma,gamma2,W_r,cE,pE,rho_E,J_sf_0,J_wf_0,Pd_0,Pd_L,Pf_L,Q_sf_0,pd_0,pf_0,pd_L,pf_L,HP_eff,LP_eff,T_eff,V_m,ERD_eff,ERD_fric,A_ERD,eta_ERD,mix_density,pw,pe,swro_beta_fix,beta_fix,mixer_ERD,version,fig,swro_KF,swro_KD,KF,KD]= Pareto_2_data(input1, input2)
+%%  Pareto_2_data(input)
+%
+%   Data for Berechung der blauen Pareto front
+%
+%   SWRO+ERD
+%
+%   option_data = 2
+%
+%   Input:
+%       input1        -   x
+%       input2        -   x
+
+%% model versions
+version=zeros(1,10);
+% version(1)=0 if co-current, 1 otherwise
+version(2)=1;  % 0 = SWRO beta fixed
+version(3)=0;  % 0 = PRO beta fixed
+version(4)=1;  % 0 = ideal SWRO
+version(5)=1;  % 0 = ideal PRO
+% ERD
+version(6)=1;
+% version(6) = 0 --> only SWRO (no ERD)
+% version(6) = 1 --> only SWRO (with ERD)
+% version(6) = 2 --> only PRO
+% version(6) = 3 --> SWRO-PRO hybrid system (with one ERD)
+% version(6) = 4 --> SWRO-PRO hybrid system (with two ERDs) 
+version(7)=1; % 1 = ICP and ECP for SWRO
+version(8)=1; % 1 = ICP and ECP for PRO
+
+%% Membrane unit properties
+H = 1e-3;           % height of the membrane [m]
+ro_water = 1000;    % mass density of water  [kg/m^3] 
+ro_salt = 2165;     % mass density of salt [kg/m^3] 
+Mw = 18;            % molecular weight of water [kg/kmol]
+Ms = 58.44;         % molecular weight of salt  [kg/kmol]
+Rw = 462;           % gas constant of water  [J/(kg K)] 
+T0 = 297;           % temperature [K] 
+eta = 1.3e-3;       % seawater viscosity [kg/(m s)]
+p_r = 1e5;          % pressure [Pa]=[kg/ms^2]
+rho_r = 1e3;        % density [kg/m^3]
+C_r = 1;            % salt concentration [%]
+
+%% SWRO
+swro_Z=1;               % width of SWRO membrane [m]
+swro_L=4;               % length of SWRO membrane [m]
+if version(6)== 2; swro_L=1; end
+swro_alpha = 5.0815e-9; % SWRO water permeablity co-efficient [s/m]
+swro_R =0.96;           % SWRO salt rejection rate 
+swro_KK = 1e-2;         % SWRO ICP mass transfer coefficient
+swro_KD = 1/swro_KK;    % SWRO ECP draw side mass transfer coefficient
+swro_KF = 1/swro_KK;    % SWRO ECP fresh side mass transfer coefficient
+swro_x_r=swro_L^2;      % x_r=swro_L^2 since x = linspace(0,1,n) (if x = linspace(0,swro_L,n) then x_r=swro_L;)     
+swro_b1 = H/swro_x_r;                           % H/swro_L ratio
+swro_b2 = swro_Z/swro_x_r;                      % Z/swro_L ratio
+J_r = sqrt(H^3/swro_x_r*p_r*rho_r);             % flux [kg/s^2]
+swro_gamma= swro_x_r * p_r *  swro_alpha /J_r;  % SWRO scaling factor - mass balance
+swro_gamma2= J_r^2./(swro_x_r^2 * p_r * rho_r); % SWRO scaling factor - momentum balance
+swro_W_r=J_r*p_r/rho_r;                         % net work [W/m^2]
+sigma=0.999 ;                                   % Rejection coefficient
+swro_beta_fix=4.43e-4/J_r*swro_x_r;             % value for fixed SWRO beta [kg/sm^2]
+ 
+%% PRO
+Z=1;                    % width of the PRO membrane [m]
+L=1.5;                  % length of the PRO membrane [m]      
+alpha = 5.47e-9;        % water permeablity co-efficient [s/m]
+R =0.94;                % salt rejection rate [1]
+KK=7.13e2;              % mass transfer coefficient [sm^2/kg]
+KD=1/KK;                % PRO ECP draw side mass transfer coefficient [sm^2/kg]
+KF=1/KK;                % PRO ECP fresh side mass transfer coefficient [sm^2/kg]
+x_r=L^2;                % x_r=L^2 since x = linspace(0,1,n) (if x = linspace(0,L,n) then x_r=L;)              
+b1 = H/x_r;             % H/L ratio
+b2 = Z/x_r;                           % Z/L ratio
+Q_r = sqrt(H^3/x_r*p_r*rho_r);        % flux [kg/s^2]
+gamma= x_r * p_r * alpha /Q_r ;       % PRO scaling factor - mass balance
+gamma2= Q_r^2./(x_r^2 * p_r * rho_r); % PRO scaling factor - momentum balance
+W_r=Q_r*p_r/rho_r;                    % net work [W/m^2]
+beta_fix =1.71e-4/Q_r*x_r;            % value for fixed PRO beta [kg/sm^2]
+
+%% Sea Water
+cE= 35/983/C_r;                            % salt concentration in seawater
+pE= 1e5/p_r;                               % external pressure
+rho_E=(cE + 1)./(ro_water*cE/ro_salt + 1); % density of incomming seawater
+
+%% SWRO operating conditions
+J_sf_0 = 0;              % salt flux in fresh side at 0 
+J_wf_0 = 0;              % water flux in fresh side at L               
+Pf_L = pE;               % pressure fresh side at L
+
+Pd_0 =  input1(1);       % pressure draw side at 0
+Pd_L =  input1(2);       % pressure draw side at L (not need in the hybrid system)
+
+%% PRO operation conditions 
+Q_sf_0 = 0;              % salt flux in fresh side at 0
+pf_L = pE;               % pressure of fresh side at L
+
+pd_0 = 12e5/p_r;         % pressure draw side at 0
+pd_L = 14.6e5/p_r;       % pressure of fresh side at 0
+pf_0 = 1.3e5/p_r;        % pressure draw side at L 
+
+%% ERD/Turbine/Pump parameters
+T_eff  = .95;               % turbine efficiency
+HP_eff = .9;            	% high pressure pump efficiency
+LP_eff = .95;               % low pressure pump efficiency
+V_m=0.052;                  % Volumetric mixing
+ERD_eff=.96;                % ERD unit pressure efficiency
+ERD_fric=5e-04;             % ERD friction coefficient
+A_ERD=H*swro_Z;             % cross sectional area of ERD inflows/outflows
+eta_ERD=0.01;               % leak of high pressure brine
+mix_density=997/rho_r;  	% density of mixture in ERD
+pw=2.4;                     % water price [$/m^3]
+pe=0.14;                    % electricity price [$/kWh]
+mixer_ERD=1;                % PRO Draw outlet mixer adjustment (only if 2nd ERDs) (mixer_ERD=1 --> all flow to ERD2 no turbine needed)    
+
+%% display figures
+%fig=[1,1,1,0]; % f(i)=1 --> figure i will be displayed
+fig=[0,0,0,0]; % f(i)=1 --> figure i will be displayed
+
+%% model specific changes:
+% co-current
+if pd_0 > pd_L; version(1)=0; else; version(1)=1; end            
+% no SWRO needed (--> use trivial data)
+if version(6) == 2; version(4)=0; Pd_0=40; Pd_L=35; end
+
+end
\ No newline at end of file
diff --git a/Sim_0a_data.m b/Sim_0a_data.m
index bc0f079..4ef6510 100644
--- a/Sim_0a_data.m
+++ b/Sim_0a_data.m
@@ -108,13 +108,13 @@ ERD_fric=5e-04;             % ERD friction coefficient
 A_ERD=H*swro_Z;             % cross sectional area of ERD inflows/outflows
 eta_ERD=0.01;               % leak of high pressure brine
 mix_density=997/rho_r;  	% density of mixture in ERD
-pw=2;                       % water price [$/m^3]
-pe=0.5;                     % electricity price [$/kWh]
+pw=2.4;                     % water price [$/m^3]
+pe=0.14;                    % electricity price [$/kWh]
 mixer_ERD=1;                % PRO Draw outlet mixer adjustment (only if 2nd ERDs) (mixer_ERD=1 --> all flow to ERD2 no turbine needed)    
 
 %% display figures
-fig=[1,0,0,0,0,0,0,0,1,0,0]; % f(i)=1 --> figure i will be displayed
-%fig=[0,0,0,0,0,0,0,0,0,0,0]; % f(i)=1 --> figure i will be displayed
+fig=[1,1,1,0]; % f(i)=1 --> figure i will be displayed
+%fig=[0,0,0,0]; % f(i)=1 --> figure i will be displayed
 
 %% model specific changes:
 % co-current
diff --git a/Sim_0b_data.m b/Sim_0b_data.m
index 2cd106e..900ff9e 100644
--- a/Sim_0b_data.m
+++ b/Sim_0b_data.m
@@ -108,13 +108,13 @@ ERD_fric=5e-04;             % ERD friction coefficient
 A_ERD=H*swro_Z;             % cross sectional area of ERD inflows/outflows
 eta_ERD=0.01;               % leak of high pressure brine
 mix_density=997/rho_r;  	% density of mixture in ERD
-pw=2;                       % water price [$/m^3]
-pe=0.5;                     % electricity price [$/kWh]
+pw=2.4;                     % water price [$/m^3]
+pe=0.14;                    % electricity price [$/kWh]
 mixer_ERD=1;                % PRO Draw outlet mixer adjustment (only if 2nd ERDs) (mixer_ERD=1 --> all flow to ERD2 no turbine needed)    
 
 %% display figures
-%fig=[1,0,0,0,0,0,0,0,1,0,0]; % f(i)=1 --> figure i will be displayed
-fig=[0,0,0,0,0,0,0,0,0,0,0]; % f(i)=1 --> figure i will be displayed
+%fig=[1,1,1,0]; % f(i)=1 --> figure i will be displayed
+fig=[0,0,0,0]; % f(i)=1 --> figure i will be displayed
 
 %% model specific changes:
 % co-current
diff --git a/Skript.asv b/Skript.asv
index e92418f..7304770 100644
--- a/Skript.asv
+++ b/Skript.asv
@@ -11,7 +11,7 @@
 % Create figures form chapter: Numerical simulations
 % 
 close all;
-[Output1, Output2, time]=fun_1(1,0,'sol',1e4,1e-4),
+[Output1, Output2, time]=fun_1(1,0,'sol',1e4,1e-6),
 %
 %figure(1); set(gcf,'color','w'); f = gcf; exportgraphics(f,'Figure_6.png');
 %figure(2); set(gcf,'color','w'); f = gcf; exportgraphics(f,'Figure_7.png');
@@ -20,24 +20,78 @@ close all;
 
 %% Average simulation time:
 %
-% IMPORTANT: manually comment out line: 116<->117 in "SIM_0_data"
+% calculate average computation time of the simulation above
 %
 T=zeros(1,100); 
 for i=1:length(T)
-[Output1, Output2, time]=fun_1(1,0,'sol',1e4,1e-4);
+[~, ~, time]=fun_1(1,0.1,'sol',1e4,1e-6);
 T(i)=time;
 end
 mean(T) %= 0.9600
 %
 
-%% Optimisation - case2 (with epsilon constraints)
+%% test 
+fun_1([70 50],option_data,'sol',option_mesh,option_BVP),
+
+%% Optimisation - case1 (epsilon constraint method)
+%
+% Create Paretofront for SWRO+ERD 
+%
+load("DATA_case1.mat")
+% 
+% load("DATA_case1.mat")
+% save DATA_case1.mat X1 Y1 X1_P Y1_P time1 time1_P
+% 
+% X1  -  set of operating conditions for epsilon constraint method (2x40)
+% Y1  -  system output using X2 (5x40)
+%
+% X1_P  -  set of operating conditions from paretosearch (2x200)
+% Y1_P  -  system output using X2_P (4x200)
+%
+% time1    - time for epsilon constraint method (40points)
+% time1_P  - time for paretosearch solver (200 points)
+%
+startTime=datetime("now");
+A= [-1 1]; b= -.1; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70];
+option_mesh = 1e4; option_BVP = 1e-6; option_data = 1;
+foptions = optimoptions('fmincon','Display','off','Algorithm','interior-point', 'StepTolerance', 1e-12, 'OptimalityTolerance',1e-3, 'MaxFunEvals',100);
+%foptions = optimoptions('fmincon','Display','iter-detailed','Algorithm','interior-point', 'StepTolerance', 1e-12, 'OptimalityTolerance',1e-3, 'MaxFunEvals',100);
+rng default
+%
+FWmin=linspace(0.0065, 0.3861, 40); 
+X0=[linspace(60,69.99,40);linspace(55,50,40)];
+%
+parfor i=1:40
+    i,
+    epsilon=1;
+    minFW=-FWmin(i);
+    x0=X0(:,i);
+    while epsilon > .9
+        [x_neu, minSEC] = fmincon(@(x)fun_1(x,option_data,'SEC',option_mesh,option_BVP),x0   ,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'FW' ,option_data,option_mesh,option_BVP,-minFW ),foptions);
+        [x0, minFW]     = fmincon(@(x)fun_1(x,option_data,'FW' ,option_mesh,option_BVP),x_neu,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'SEC',option_data,option_mesh,option_BVP,-minSEC),foptions);
+        epsilon = max(norm(x0-x_neu)),
+    end
+    X1(:,i)=x0; %x_neu;
+    Y1(:,i)=fun_1([X1(:,i)],option_data,'sol',option_mesh,option_BVP);
+    clc,
+end
+endTime=datetime("now");
+time1 = endTime - startTime;
+save DATA_case1.mat X1 Y1 X1_P Y1_P time1 time1_P
+% scatter
+f=figure(1); f.Position = [1200 500 800 500];
+scatter(Y1(1,:),Y1(2,:),'MarkerEdgeColor',[0 0 0],'MarkerFaceColor','r'); hold on
+grid on;lgd=legend('Single SWRO unit','Location', 'Northeast');ylabel('FW [m^3/h]','FontSize',16);xlabel('SEC_{net} [kWh/m^3]','FontSize',16); %,'SWRO with ERD', 'Hybrid system', 
+xlim([-5.501 -.5]); ylim([0 0.48]);view(2);
+
+%% Optimisation - case2 (epsilon constraint method) -- FIX Y dimension to 5!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 %
 % Create Paretofront for SWRO+ERD 
 %
 load("DATA_case2.mat")
 % 
 % load("DATA_case2.mat")
-% save DATA_case2.mat X_2
+% save DATA_case2.mat X2 Y2 X2_P Y2_P time2 time2_P
 % 
 % X2  -  set of operating conditions for epsilon constraint method (2x40)
 % Y2  -  system output using X2 (4x40)
@@ -48,31 +102,66 @@ load("DATA_case2.mat")
 % time2    - time for epsilon constraint method (40points)
 % time2_P  - time for paretosearch solver (200 points)
 %
+startTime=datetime("now");
+A= [-1 1]; b= -.1; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70];
+option_mesh = 1e4; option_BVP = 1e-6; option_data = 1;
+foptions = optimoptions('fmincon','Display','off','Algorithm','interior-point', 'StepTolerance', 1e-12, 'OptimalityTolerance',1e-3, 'MaxFunEvals',100);
+rng default
+%
+FWmin=linspace(0.0017, 0.37, 40); 
+X0=[linspace(37,70,40);linspace(36,50,40)];
+%
+parfor i=1:40
+    i,
+    epsilon=1;
+    minFW=-FWmin(i);
+    x0=X0(:,i);
+    while epsilon > .9
+        [x_neu, minSEC] = fmincon(@(x)fun_1(x,option_data,'SEC',option_mesh,option_BVP),x0   ,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'FW' ,option_data,option_mesh,option_BVP,-minFW ),foptions);
+        [x0, minFW]     = fmincon(@(x)fun_1(x,option_data,'FW' ,option_mesh,option_BVP),x_neu,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'SEC',option_data,option_mesh,option_BVP,-minSEC),foptions);
+        epsilon = max(norm(x0-x_neu)),
+    end
+    X1(:,i)=x0; %x_neu;
+    Y1(:,i)=fun_1([X1(:,i)],option_data,'sol',option_mesh,option_BVP);
+    clc,
+end
+endTime=datetime("now");
+time2 = endTime - startTime;
+save DATA_case2.mat X2 Y2 X2_P Y2_P time2 time2_P
+%scatter
+scatter(Y2(1,:),Y2(2,:),'MarkerEdgeColor',[0 0 0],'MarkerFaceColor','b'); hold on
+
+
+
+
+
+
+
+
+
+
+
+
 
-X2=zeros(2,40);Y2=zeros(4,40);X2_P=zeros(2,200);
 
 
 
-load("DATA_T1.mat")
-c1=datetime("now");
-x0=[69.9 51.5];
-FWmin=linspace(0.3770, 0.0015, 40);
-for i=1:40
-clc,
-i, 
-close all; A= [-1 1]; b= -.1; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70];
-foptions = optimoptions('fmincon','Display','iter-detailed','Algorithm','interior-point', 'StepTolerance', 1e-16, 'OptimalityTolerance',1e-4, 'MaxFunEvals',100);
-option_mesh = 1e4; option_BVP = 1e-4; option_data = 1;
-[x_opt, f_opt]=fmincon(@(x)fun_1(x,1,'SEC',option_mesh,option_BVP),x0,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'FW',1,option_mesh,option_BVP, FWmin(i)),foptions);
-end
-c2=datetime("now"); %von 2 bis 39
-save DATA_T1.mat c1 c2 % 00:02:30
 
 %% Optimisation - case2 (with paretosearch)
+load("DATA_case2.mat")
 %
-% Create Paretofront for SWRO+ERD 
+startTime=datetime("now");
+A= [-1 1]; b= -.1; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70];
+option_mesh = 1e4; option_BVP = 1e-4; option_data = 2;
 %
-
+rng default
+% take 200 random data points
+% run paretosearch 
+%
+% endTime
+% runTime ...
+% scater plot togerther with other 40 points
+% define all variables as floating points and see that the
 
 
 %% Optimisation-2 (red)
@@ -123,6 +212,10 @@ scatter(SEC_Pareto_2,FW_Pareto_2)
 
 %% red - Paretosearch
 %
+%  UseCompletePoll =true --> paretosearch in parallel
+%  UseParallel =true --> paralell paretosearch
+%
+%
 close all; clear all; clc
 rng default % For reproducibility
 %
diff --git a/Skript.m b/Skript.m
index 2a11bef..b820493 100644
--- a/Skript.m
+++ b/Skript.m
@@ -11,7 +11,7 @@
 % Create figures form chapter: Numerical simulations
 % 
 close all;
-[Output1, Output2, time]=fun_1(1,0,'sol',1e4,1e-4),
+[Output1, Output2, time]=fun_1(1,0,'sol',1e4,1e-6),
 %
 %figure(1); set(gcf,'color','w'); f = gcf; exportgraphics(f,'Figure_6.png');
 %figure(2); set(gcf,'color','w'); f = gcf; exportgraphics(f,'Figure_7.png');
@@ -20,17 +20,66 @@ close all;
 
 %% Average simulation time:
 %
-% calculate average computation time of previous simulation
+% calculate average computation time of the simulation above
 %
 T=zeros(1,100); 
 for i=1:length(T)
-[Output1, Output2, time]=fun_1(1,0.1,'sol',1e4,1e-4);
+[~, ~, time]=fun_1(1,0.1,'sol',1e4,1e-6);
 T(i)=time;
 end
 mean(T) %= 0.9600
 %
 
-%% Optimisation - case2 (with epsilon constraints)
+%% Optimisation - case1 (epsilon constraint method)
+%
+% Create Paretofront for SWRO+ERD 
+%
+load("DATA_case1.mat")
+% 
+% load("DATA_case1.mat")
+% save DATA_case1.mat X1 Y1 X1_P Y1_P time1 time1_P
+% 
+% X1  -  set of operating conditions for epsilon constraint method (2x40)
+% Y1  -  system output using X2 (5x40)
+%
+% X1_P  -  set of operating conditions from paretosearch (2x200)
+% Y1_P  -  system output using X2_P (4x200)
+%
+% time1    - time for epsilon constraint method (40points)
+% time1_P  - time for paretosearch solver (200 points)
+%
+startTime=datetime("now");
+A= [-1 1]; b= -.1; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70];
+option_mesh = 1e4; option_BVP = 1e-6; option_data = 1;
+foptions = optimoptions('fmincon','Display','off','Algorithm','interior-point', 'StepTolerance', 1e-12, 'OptimalityTolerance',1e-4, 'MaxFunEvals',100);
+%foptions = optimoptions('fmincon','Display','iter-detailed','Algorithm','interior-point', 'StepTolerance', 1e-12, 'OptimalityTolerance',1e-3, 'MaxFunEvals',100);
+rng default
+%
+FWmin=linspace(0.0065, 0.3861, 40); 
+X0=[linspace(60,69.99,40);linspace(55,50,40)];
+%
+parfor i=1:40
+    epsilon=1;
+    minFW=-FWmin(i);
+    x0=X0(:,i);
+    while epsilon > .9
+        [x_neu, minSEC] = fmincon(@(x)fun_1(x,option_data,'SEC',option_mesh,option_BVP),x0   ,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'FW' ,option_data,option_mesh,option_BVP,-minFW ),foptions);
+        [x0, minFW]     = fmincon(@(x)fun_1(x,option_data,'FW' ,option_mesh,option_BVP),x_neu,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'SEC',option_data,option_mesh,option_BVP,-minSEC),foptions);
+        epsilon = max(norm(x0-x_neu)),
+    end
+    X1(:,i)=x0; %x_neu;
+    Y1(:,i)=fun_1([X1(:,i)],option_data,'sol',option_mesh,option_BVP);
+end
+endTime=datetime("now");
+time1 = endTime - startTime;
+save DATA_case1.mat X1 Y1 X1_P Y1_P time1 time1_P
+%scatter
+f=figure(1); f.Position = [1200 500 800 500];
+scatter(Y1(1,:),Y1(2,:),'MarkerEdgeColor',[0 0 0],'MarkerFaceColor','r'); hold on
+grid on; xlim([-5.501 -.5]); ylim([0 0.48]);view(2);
+legend('Single SWRO unit','SWRO with ERD','Location', 'Northeast');ylabel('FW [m^3/h]','FontSize',16);xlabel('SEC_{net} [kWh/m^3]','FontSize',16);
+
+% Optimisation - case2 (epsilon constraint method)
 %
 % Create Paretofront for SWRO+ERD 
 %
@@ -48,27 +97,78 @@ load("DATA_case2.mat")
 % time2    - time for epsilon constraint method (40points)
 % time2_P  - time for paretosearch solver (200 points)
 %
-
-load("DATA_T1.mat")
-c1=datetime("now");
-x0=[69.9 51.5];
-FWmin=linspace(0.3770, 0.0015, 40);
-for i=1:40
-clc,
-i, 
-close all; A= [-1 1]; b= -.1; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70];
-foptions = optimoptions('fmincon','Display','iter-detailed','Algorithm','interior-point', 'StepTolerance', 1e-16, 'OptimalityTolerance',1e-4, 'MaxFunEvals',100);
-option_mesh = 1e4; option_BVP = 1e-4; option_data = 1;
-[x_opt, f_opt]=fmincon(@(x)fun_1(x,1,'SEC',option_mesh,option_BVP),x0,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'FW',1,option_mesh,option_BVP, FWmin(i)),foptions);
+startTime=datetime("now");
+A= [-1 1]; b= -.1; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70];
+option_mesh = 1e4; option_BVP = 1e-6; option_data = 1;
+foptions = optimoptions('fmincon','Display','off','Algorithm','interior-point', 'StepTolerance', 1e-12, 'OptimalityTolerance',1e-4, 'MaxFunEvals',100);
+rng default
+%
+FWmin=linspace(0.0015, 0.3769, 40); 
+X0=[linspace(37,70,40);linspace(36,50,40)];
+%
+parfor i=1:40
+    epsilon=1;
+    minFW=-FWmin(i);
+    x0=X0(:,i);
+    while epsilon > .9
+        [x_neu, minSEC] = fmincon(@(x)fun_1(x,option_data,'SEC',option_mesh,option_BVP),x0   ,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'FW' ,option_data,option_mesh,option_BVP,-minFW ),foptions);
+        [x0, minFW]     = fmincon(@(x)fun_1(x,option_data,'FW' ,option_mesh,option_BVP),x_neu,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'SEC',option_data,option_mesh,option_BVP,-minSEC),foptions);
+        epsilon = max(norm(x0-x_neu)),
+    end
+    X2(:,i)=x0; %x_neu;
+    Y2(:,i)=fun_1([X1(:,i)],option_data,'sol',option_mesh,option_BVP);
 end
-c2=datetime("now"); %von 2 bis 39
-save DATA_T1.mat c1 c2 % 00:02:30
+endTime=datetime("now");
+time2 = endTime - startTime;
+save DATA_case2.mat X2 Y2 X2_P Y2_P time2 time2_P
+%scatter
+scatter(Y2(1,:),Y2(2,:),'MarkerEdgeColor',[0 0 0],'MarkerFaceColor','b'); hold on
+grid on; xlim([-5.501 -.5]); ylim([0 0.48]);view(2);
+legend('Single SWRO unit','SWRO with ERD','Location', 'Northeast');ylabel('FW [m^3/h]','FontSize',16);xlabel('SEC_{net} [kWh/m^3]','FontSize',16);
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 
 %% Optimisation - case2 (with paretosearch)
+load("DATA_case2.mat")
 %
-% Create Paretofront for SWRO+ERD 
+startTime=datetime("now");
+A= [-1 1]; b= -.1; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70];
+option_mesh = 1e4; option_BVP = 1e-4; option_data = 2;
 %
-
+rng default
+% take 200 random data points
+% run paretosearch 
+%
+% endTime
+% runTime ...
+% scater plot togerther with other 40 points
+% define all variables as floating points and see that the
 
 
 %% Optimisation-2 (red)
@@ -119,6 +219,10 @@ scatter(SEC_Pareto_2,FW_Pareto_2)
 
 %% red - Paretosearch
 %
+%  UseCompletePoll =true --> paretosearch in parallel
+%  UseParallel =true --> paralell paretosearch
+%
+%
 close all; clear all; clc
 rng default % For reproducibility
 %
diff --git a/fun_1.asv b/fun_1.asv
index 15e4f1a..9284af8 100644
--- a/fun_1.asv
+++ b/fun_1.asv
@@ -9,7 +9,7 @@ function [output1, output2, output3] = fun_1(input,option_data,obj,option_mesh,o
 %   Input:
 %       input         -   Input for the Data function
 %       option_data   -   Selects Data function used: data(input)
-%       obj           -   'SEC', 'FW', 'Rev', 'PD_net' , 'SE_net', 'approx'  or 'sol' ;
+%       obj           -   'SEC', 'FW', 'Rev', 'Pareto', 'PD_net', 'SE_net', 'SE_f', 'sol';
 %       option_mesh   -   NMax of bvp5c
 %       option_BVP    -   RelTol of bvp5c
 %
@@ -26,12 +26,17 @@ function [output1, output2, output3] = fun_1(input,option_data,obj,option_mesh,o
 %       tol=1e-6;  
 %       [Sim_01_output1, Sim_01_output2]=fun(x0,D,obj,mesh,tol),
 %
-tic;
+if nargout>2
+    tic;
+end
+
 %% Read data
 if option_data == 0; DATA = @(x)Sim_0a_data(input); end
 if option_data == 0.1; DATA = @(x)Sim_0b_data(input); end
-
-
+%
+if option_data==1; DATA = @(x)Pareto_1_data(input); end
+if option_data==2; DATA = @(x)Pareto_2_data(input); end
+%
 
 % Simulations
 if option_data == 0.11; DATA = @(x)Sim_11_data(input); end
@@ -50,13 +55,9 @@ if option_data == 0.53; DATA = @(x)Sim_PRO_03_data(input); end
 % Rough vs. Detailled
 if option_data == 0.91; DATA = @(x)Sim_rough_approx_data(input); end
 % Optimizations
-if option_data==1; DATA = @(x)Opt_1_data(input); end
-if option_data==2; DATA = @(x)Opt_2_data(input); end
 if option_data==3; DATA = @(x)Opt_3_data(input); end
 if option_data==4; DATA = @(x)Opt_4_data(input); end
 if option_data==5; DATA = @(x)Opt_L_data(input); end
-% for 3d plot
-if option_data==11; DATA = @(x)plot3d_data(input); end
 
 [H,Z,swro_Z,ro_water,ro_salt,Mw,Ms,Rw,T0,eta,sigma,p_r,rho_r,C_r,swro_L,swro_alpha,swro_R,swro_KK,swro_x_r,swro_b1,swro_b2,J_r,swro_gamma,swro_gamma2,swro_W_r,L,alpha,R,KK,x_r,b1,b2,Q_r,gamma,gamma2,W_r,cE,pE,rho_E,J_sf_0,J_wf_0,Pd_0,Pd_L,Pf_L,Q_sf_0,pd_0,pf_0,pd_L,pf_L,HP_eff,LP_eff,T_eff,V_m,ERD_eff,ERD_fric,A_ERD,eta_ERD,mix_density,pw,pe,swro_beta_fix,beta_fix,mixer_ERD,version,fig,swro_KF,swro_KD, KF, KD] ...
     = DATA(1);
@@ -114,8 +115,8 @@ if version(4) == 0 ;swro_beta=zeros(1,n); end
 % J_w,in and J_s,in
 J_cross = ((Y(:,5) - Y(:,6)) - sigma.*(swro_p_osm_d - swro_p_osm_f))./(1 + p_r*swro_alpha*swro_KK*sigma.*(swro_p_osm_d - swro_p_osm_f));
 if version(4) == 0; J_cross = ((Y(:,5) - Y(:,6)) - sigma*(swro_p_osm_d - swro_p_osm_f));  end 
-J_sin=swro_beta.*swro_del_c;
-if version(7)==1
+J_sin=swro_beta'.*swro_del_c;
+if version(7)==1 && version(4)>0
 J_cross = ((Y(:,5)-Y(:,6))-sigma.*(swro_p_osm_d - swro_p_osm_f) + swro_beta.*(Y(:,5)-Y(:,6)).*(swro_KK +1/swro_KD + 1/swro_KF)) ./ (1 + swro_alpha.*swro_beta.*(1/swro_KD + 1/swro_KF + swro_KK) + p_r*swro_alpha*sigma.*(-swro_p_osm_d./swro_KD -swro_p_osm_f.*(swro_KK + 1/swro_KF) ));
 J_sin = swro_beta.*( swro_del_c-J_cross.*( Y(:,1)./(Y(:,1) + 1)./swro_KD +  Y(:,3)./(Y(:,3) + Y(:,4)).*(swro_KK+1/swro_KF)))./(1+swro_alpha.*swro_beta.*(swro_KK+1/swro_KF+1/swro_KD));
 end
@@ -144,7 +145,7 @@ if version(5) == 0 ; beta = zeros(1,n); end
 Q_cross = (p_osm_d - p_osm_f - (Y(:,5+6) - Y(:,6+6)) - (Y(:,5+6) - Y(:,6+6)).*beta'.*KK.*alpha.*p_r)./((ones(1,n) + KK*alpha*p_r.*(beta +p_osm_f'))');
 if version(5) == 0; Q_cross = ((p_osm_d - p_osm_f) - (Y(:,5+6) - Y(:,6+6))) ; end
 Q_sin = beta'.*del_c;
-if version(8)==1 
+if version(8)==1 && version(4)>0
 Q_cross = ((p_osm_d - p_osm_f) - (Y(:,5+6) - Y(:,6+6)) - (Y(:,5+6) - Y(:,6+6)).*beta'.*(KK+1./KF+1./KD).*p_r.*alpha)./(ones(1,n) + alpha.*beta.*(KK+1./KF+1./KD) + alpha.*p_r.*(p_osm_f'.*(KK+1./KF) + p_osm_d'./KD))';
 Q_sin = beta'.*( del_c-Q_cross.*( Y(:,1+6)./(Y(:,1+6) + 1)./KD +  Y(:,3+6)./(Y(:,3+6) + Y(:,4+6)).*(KK+1/KF)))./(1+alpha.*beta');
 end
@@ -164,13 +165,14 @@ PRO_Recovery = (1- Q_f(end)./Q_f(1))*100;  % PRO recovery rate [%]
     if SWRO_Recovery < 0 
         fprintf(2,' \nERROR: Waring: SWRO recovery is negative! \n');
     end
-    if PRO_Recovery > 100 
-        fprintf(2,' \nERROR: Waring: PRO recovery is greater then 100 %% \n');
-    end
-    if PRO_Recovery < 0 
-        fprintf(2,' \nERROR: Waring: PRO recovery is negative! \n');
+    if version(6)>1
+        if PRO_Recovery > 100 
+            fprintf(2,' \nERROR: Waring: PRO recovery is greater then 100 %% \n');
+        end
+        if PRO_Recovery < 0 
+            fprintf(2,' \nERROR: Waring: PRO recovery is negative! \n');
+        end
     end
-
 %% Output of System
 J_E = (C_d(1)*J_wd(1)+J_wd(1))/2; J_ERD = J_E; J_wE = J_E/(cE+1);
 W_p1 = 1/HP_eff * (P_d(1)-pE)*(J_E *swro_Z)/rho_E; W_p1 = W_p1*swro_W_r;
@@ -192,11 +194,15 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r;
         f_2= (J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2); 
         f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
         pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_exit*Q_exit*swro_Z/rho_exit ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;
+        pMax = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_exit*Q_exit*swro_Z/rho_exit ) +  pE*J_E*swro_Z/rho_E )/J_ERD/swro_Z;
+        if pERD > pMax; pERD2=pMax; end
+        if pERD < 0; pERD2=pMax; end  
+
         W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD; W_p3=W_p3*swro_W_r; W_p4=0; W_t=0;
             if W_p3 < 0 
                 fprintf(2,' \nERROR: Waring: Pump 3 generates generates power! \n');
+                sol.stats.maxerr =1000;
             end
-
     end
     %% Version(6)=2 --> only PRO
     if version(6)==2
@@ -218,6 +224,8 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r;
                 f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2;  
                 f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
                 pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;
+                pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pE*J_E*swro_Z/rho_E)/J_ERD/swro_Z;
+                pERD = min(max(1,pERD),pMax);
                 W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD;
                 W_t = T_eff * (p_d(end)-pE)*(Q_d(end)*Z)/local_ro_d(end); 
             end
@@ -226,6 +234,8 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r;
                 f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2  ;
                 f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
                 pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;     
+                pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pE*J_E*swro_Z/rho_E )/J_ERD/swro_Z;       
+                pERD = min(max(1,pERD),pMax);
                 W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD;
                 W_t = T_eff * (p_d(1)-pE)*(abs(Q_d(1))*Z)/local_ro_d(1);
             end
@@ -251,52 +261,6 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r;
             Q_exit = mixer_ERD*abs(Q_d(1))*(1-eta_ERD);
             c_exit = (c_d(1)*mixer_ERD*abs(Q_wd(1))*(1 - eta_ERD)+cE*J_wE_sea-C_ERD2*2*J_E)/(mixer_ERD*abs(Q_wd(1))*(1 - eta_ERD)+J_wE_sea-2*J_E);
         end
-
-
-%% Removed comments
-
-%Y(:,4+6)=abs(Y(:,4+6)) % ??why is this here??
-
-% -3.9191    0.3247    0.0131   32.8727    7.1201       NaN       NaN
-%
-% 1.0e+03 *
-% -1.2742   -0.9744         0   -0.4873   -0.0010    0.1885
-
-% In ideal cases take mean of some quantities. 
-% In ideal cases salt will fluctuate around 1e-16. We take mean here for 
-% 
-% This is just for figures. (), otherwise one could use ylim ) 
-%if version(4)==0 ; J_f=J_wf; J_sf=zeros(1,n); J_sd=mean(J_sd)*ones(1,n); C_f=zeros(1,n); end
-%if version(5)==0 ; Q_f=Q_wf; Q_sf=zeros(1,n); Q_sd=mean(Q_sd)*ones(1,n); c_f=zeros(1,n); end
-
-        %f_ERD=0; pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_exit*Q_exit*swro_Z/rho_exit ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;
-        %if pERD > pMax; pERD=pMax; end
-        %if pERD < 0; pERD=pMax; end  
-        %display(['pERD:              ',num2str(pERD)]) % pERD use in the model
-        %display(['pERD_theoretical:  ',num2str(pMax)]) % theoretical pERD with no leak, friction perm and 100% efficiency
-
-        %f_ERD=0;pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;
-        %if pERD > pMax; pERD=pMax; end
-        %if pERD < 0; pERD=pMax; end  
-
-        %f_ERD=0;pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;       
-        %display(['pERD:              ',num2str(pERD)]) % pERD use in the model
-        %display(['pERD_theoretical:  ',num2str(pMax)]) % theoretical pERD with no leak, friction perm and 100% efficiency
-        %if pERD > pMax; pERD=pMax; end
-        %if pERD < 0; pERD=pMax; end     
-
-        %J_wERD = J_E/(C_ERD+1);
-
-
-
-
-
-
-
-
-
-
-
         rho_exit= (ro_salt*(c_exit+1))/(c_exit*ro_water+ro_salt); p_exit=pE;
         %calculation of Pressure P_ERD2
             if version(1) == 0 %co-current
@@ -304,19 +268,13 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r;
                 f_2= J_E*swro_Z/rho_ERD2*mix_density*(qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end)/A_ERD)^2;  
                 f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
                 pERD2 = rho_ERD2*(ERD_eff*( p_d(end)*qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z;
-                f_ERD=0;pMax = rho_ERD2*(ERD_eff*( p_d(end)*qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z;
-                if pERD2 > pMax; pERD2=pMax; end
-                %if pERD2 < 0; pERD2=pMax; end  
-                
             end
             if version (1) == 1 %counter-current 
                 f_1= qq*Q_exit*Z/rho_exit*mix_density*(J_E_sea*swro_Z/rho_E/A_ERD)^2;  
                 f_2= J_E*swro_Z/rho_ERD2*mix_density*(qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1)/A_ERD)^2;  
                 f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
                 pERD2 = rho_ERD2*(ERD_eff*( p_d(1)*qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z;            
-                f_ERD=0;pMax = rho_ERD2*(ERD_eff*( p_d(1)*qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z;            
-                if pERD2 > pMax; pERD2=pMax; end
-                %if pERD2 < 0; pERD2=pMax; end  
+
             end     
             %calculation of Pressure P_ERD
             if version(1) == 0 %co-current
@@ -324,9 +282,6 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r;
                 f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2;  
                 f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
                 pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;
-                f_ERD=0; pMax = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;
-                if pERD > pMax; pERD=pMax; end
-                %if pERD < 0; pERD=pMax; end  
                 W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD;
                 W_t = T_eff * (p_d(end)-pE)*((1-mixer_ERD)*Q_d(end)*Z)/local_ro_d(end); 
             end
@@ -334,60 +289,62 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r;
                 f_1= abs(qq*Q_d(end))*Z/local_ro_d(end)*mix_density*(J_E*swro_Z/rho_ERD2/A_ERD)^2  ;
                 f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2  ;
                 f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
-                pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;            
-                f_ERD=0; pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;            
-                if pERD > pMax; pERD=pMax; end 
-                %if pERD < 0; pERD=pMax; end  
+                pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;           
                 W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD;
                 W_t = T_eff * (p_d(1)-pE)*(abs((1-mixer_ERD)*Q_d(1))*Z)/local_ro_d(1);
-            end         
+            end
     % Recalculation for Pump 1
     W_p1 = 1/HP_eff * (P_d(1)-pERD2)*(J_E *swro_Z)/rho_ERD2; W_p1 = W_p1*swro_W_r;
     end
     
 %% final output:
 W_net= - W_p1 - W_p2 - W_p3 - W_p4 + W_t; % in [W]
-PD_net = W_net./(Z*L); %in [W/m^2]    
-%old SE_net
-%if version(1)==0; SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)) + Q_d(1)*Q_r*Z/(rho_r*local_ro_d(1)))/1000/3600; end %in [kWh/m^3]
-%if version(1)==1; SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)) + abs(Q_d(end))*Q_r*Z/(rho_r*local_ro_d(end)))/1000/3600; end %in [kWh/m^3]
-%new SE_tot
-if version(1)==0; SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)))/1000/3600; end %in [kWh/m^3]
-if version(1)==1; SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)))/1000/3600; end %in [kWh/m^3]
+
+if version(6)>1
+PD_net = W_net./(Z*L); %in [W/m^2]   
+    if version(1)==0 
+        SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)) + Q_d(1)*Q_r*Z/(rho_r*local_ro_d(1)))/1000/3600; %in [kWh/m^3] 
+        SE_f   = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)))/1000/3600; %in [kWh/m^3]
+    else 
+        SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)) + abs(Q_d(end))*Q_r*Z/(rho_r*local_ro_d(end)))/1000/3600; %in [kWh/m^3] 
+        SE_f   = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)))/1000/3600; %in [kWh/m^3]
+    end 
+end
 SEC_net = W_net*(swro_local_ro_f(end)*rho_r)./(J_f(end)*J_r*swro_Z)/1000/3600; %in [kWh/m^3]  
+    if SEC_net > 0 
+        fprintf(2,' \nERROR: Waring: SEC_net is positive! \n');
+    end
 Rev= pw*FW + pe*SEC_net*FW; %in [$/h]
-tt=toc;
+if nargout>2
+    output3=toc;
+end
+
 %% if BVP-solver failed
 catch
     sol.stats.maxerr =1000;
-    SEC_net=1000;FW=0;Rev=-20;SWRO_Recovery=NaN;PRO_Recovery=NaN;PD_net=NaN;SE_net=NaN;
-    W_net=NaN;W_p1=NaN;W_p2=NaN;W_p3=NaN;W_p4=NaN;W_t=NaN;
 end
+
 %% Output of objective function
 switch (obj)
-        case 'SEC'
+        case 'SEC' % for SEC_net maximization  -->  output1 = -SEC_net
             if sol.stats.maxerr > option_BVP
-            output1 = 20;
+            output1 = 20; % using NaN causes "fmincon" to stop prematurely --> alternative "patternsearch"
             else; output1 = -SEC_net;
             end 
-        case 'FW' 
+        case 'FW' % for FW maximization  -->  output1 = -FW
             if sol.stats.maxerr > option_BVP
-            output1 = 20;
+            output1 = 20; 
             else; output1 = -FW;
             end 
-        case 'Rev' 
+        case 'Rev' % for Rev maximization  -->  output1 = -Rev
             if sol.stats.maxerr > option_BVP
             output1 = 0;
             else; output1 = -Rev;
             end
-        case 'Pareto' 
+        case 'Pareto' % SEC_net and FW maximization  -->  output1 = [-SEC_net, -FW]
             if sol.stats.maxerr > option_BVP    
-            output1 = [NaN, NaN];
+            output1 = [20, 20]; %[NaN, NaN];
             else; output1 = [-SEC_net, -FW];
-                %if SWRO_Recovery < 30;output1 = [NaN, NaN];end
-                %if SWRO_Recovery > 60;output1 = [NaN, NaN];end
-                %if PRO_Recovery < 40;output1 = [NaN, NaN];end
-                %if PRO_Recovery > 95;output1 = [NaN, NaN];end
             end 
         case 'PD_net' 
             if sol.stats.maxerr > option_BVP
@@ -399,54 +356,25 @@ switch (obj)
             output1 = 0;
             else; output1 = -SE_net;
             end
-        case 'approx' 
-            if version(6) == 0; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end
-            if version(6) == 1; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end
-            if version(6) == 2; output1=[NaN, NaN, NaN, NaN, PRO_Recovery, PD_net, SE_net]; end
-            if version(6) == 3; output1=[SEC_net, FW, Rev, SWRO_Recovery, PRO_Recovery, NaN, NaN]; end
-            if version(6) == 4; output1=[SEC_net, FW, Rev, SWRO_Recovery, PRO_Recovery, NaN, NaN]; end
-            %sol.stats.nmeshpoints
-            %sol.stats.maxerr
-            %sol.stats.nODEevals
-            %sol.stats.nBCevals
-            output2=[sol.stats.nmeshpoints sol.stats.maxerr];
-    case 'exe' 
+        case 'SE_f' 
             if sol.stats.maxerr > option_BVP
-            output1 = [NaN, NaN, NaN, NaN, NaN, NaN, NaN];
-            output2=[NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,NaN, NaN, NaN, NaN, NaN, NaN,NaN, NaN, NaN, NaN, NaN, NaN];
-            else
-            if version(6) == 0; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end
-            if version(6) == 1; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end
-            if version(6) == 2; output1=[NaN, NaN, NaN, NaN, PRO_Recovery, PD_net, SE_net]; end
-            if version(6) == 3; output1=[SEC_net, FW, SWRO_Recovery, PRO_Recovery]; end
-            if version(6) == 4; output1=[SEC_net, FW, SWRO_Recovery, PRO_Recovery]; end
-            output2=[x*swro_L; 100*C_d; J_r*J_d; p_r*P_d; p_r*P_f; x*L; 100*c_d; Q_r*Q_d; Q_r*Q_f; p_r*p_d; p_r*p_f]; % 11
-            output3=[x(2:end)*swro_L; 100*C_f(2:end); J_r*J_f(2:end); J_r/swro_x_r*J_cross(2:end); J_r/swro_x_r*J_sin(2:end); x(2:end)*L; 100*c_f(2:end); Q_r/x_r*Q_cross(2:end); Q_r/x_r*Q_sin(2:end)]; %9
-            %
-            output1=output1';output2=output2';output3=output3';
-            [a,b]=size(output1);
-            output1=[1:1:b; output1];
-            [a,b]=size(output2);
-            output2=[1:1:b; output2];
-            [a,b]=size(output3);
-            output3=[1:1:b; output3];
-            end
-    case 'sol' 
+            output1 = 0;
+            else; output1 = -SE_f;
+            end    
+        case 'sol' 
             if sol.stats.maxerr > option_BVP
-            output1 = [NaN, NaN, NaN, NaN, NaN, NaN, NaN];
-            output2=[NaN, NaN, NaN, NaN, NaN, NaN];
-            output3=NaN;
+            output1 = [NaN, NaN, NaN, NaN, NaN];
+            output2 = [NaN, NaN, NaN, NaN, NaN, NaN];
+            output3 = NaN;
             else
-            if version(6) == 0; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end
-            if version(6) == 1; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end
-            if version(6) == 2; output1=[NaN, NaN, NaN, NaN, PRO_Recovery, PD_net, SE_net]; end
-            if version(6) == 3; output1=[SEC_net, FW, Rev, SWRO_Recovery, PRO_Recovery, NaN, NaN]; end
-            if version(6) == 4; output1=[SEC_net, FW, Rev, SWRO_Recovery, PRO_Recovery, NaN, NaN]; end
+            if version(6) == 0; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN]; end
+            if version(6) == 1; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN]; end
+            if version(6) == 2; output1=[PRO_Recovery, PD_net, SE_net, SE_f, NaN]; end
+            if version(6) == 3; output1=[SEC_net, FW, Rev, SWRO_Recovery, PRO_Recovery]; end
+            if version(6) == 4; output1=[SEC_net, FW, Rev, SWRO_Recovery, PRO_Recovery]; end
             output2=[W_net, -W_p1, -W_p2,  -W_p3, -W_p4, W_t];
-            output3=tt;
             end
-  
-            
+      
 %% figure 1
 if sol.stats.maxerr ==1000
 fprintf(2,' \nERROR: bvp5c could not satisfy the relative error tolerance ---> no figures could be displayed\n');
@@ -486,13 +414,16 @@ subplot(3,2,6);lc='#0072bD';rc='#77AC30';
 plot(x*L, p_d*p_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[Pa]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 5; hold on
 yyaxis right
 plot(x*L, p_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{PRO}(x)','P_f^{PRO}(x)','Location','northeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+end
+%% figure 2
+if fig(2) == 1
 f=figure(2); %
 f.Position = [199 232.3333 1.0987e+03 639.3333];
 % Permeate flows 
 subplot(2,2,1); lc='k';rc='#b81414'; J_cross2=J_cross(2:end);J_sin2=J_sin(2:end);
 plot(x2*swro_L, J_cross2*J_r/swro_x_r, 'Color', lc, 'LineWidth',lw);xlabel('x [m]','Fontsize',10);  ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0;ylim([min( J_cross2*J_r/swro_x_r) 1.1*max( J_cross2*J_r/swro_x_r)]); hold on
 yyaxis right
-plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{RO}(x)','J_{s,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([min( J_sin2*J_r/swro_x_r) 1.1*max( J_sin2*J_r/swro_x_r)]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{RO}(x)','J_{s,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([min( J_sin2*J_r/swro_x_r) 1.1*max( J_sin2*J_r/swro_x_r)+.00001]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
 % Osmotic/ Hydraulic pressure difference
 subplot(2,2,3);lc='k';rc='#b81414'; osm_diff=swro_p_osm_d(2:end)-swro_p_osm_f(2:end); a=min([min(osm_diff*p_r) min((P_d-P_f)*p_r)]); b=max( [max(osm_diff*p_r), max((P_d-P_f)*p_r)]);
 plot(x*swro_L, (P_d-P_f)*p_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ylim([a b]);ay=gca; ay.YAxis.Exponent = 6; hold on
@@ -508,87 +439,9 @@ subplot(2,2,4);lc='k';rc='#b81414';osm_diff=p_osm_d(2:end)-p_osm_f(2:end); %a=mi
 plot(x*L, (p_d-p_f)*p_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ay=gca;ay.YAxis.Exponent = 6; hold on; %ylim([a b]); 
 yyaxis right; 
 plot(x2*L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{PRO}(x)','\Delta \pi^{PRO}(x)','Location','southeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; %ylim([a b]);
-
-
-
-end
-%% figure 2
-if fig(2) == 1
-figure(2); % pressure+concentration in both PRO channels
-subplot(2,2,1);plot(x*L, p_d*p_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('p_d(x)','Location','northeast'); xlim([0 L]);ay=gca; ay.YAxis.Exponent = 5;
-subplot(2,2,2);plot(x*L, p_f*p_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('p_f(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 5;
-subplot(2,2,3);plot(x*L, 100*c_d*C_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('c_d(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
-x2=x(2:end);c_f2=100.*c_f(2:end);
-subplot(2,2,4);plot(x2*L, c_f2*C_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('c_f(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
 end
 %% figure 3
 if fig(3) == 1
-figure(3); % Q_s,in and Q_w,in
-x2=x(2:end);Q_cross2=Q_cross(2:end);Q_sin2=Q_sin(2:end);
-subplot(1,2,1);plot(x2*L, Q_cross2*Q_r/x_r, 'b', 'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('Q_{w,in}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(1,2,2);plot(x2*L, Q_sin2*Q_r/x_r, 'g', 'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('Q_{s,in}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
-end
-%% figure 4
-lw=1.5; %Linewidth for all figures
-if fig(4) == 1
-figure(4); % Beta of PRO
-    if version(2)==0
-    plot(x*swro_L, swro_beta*J_r/swro_x_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('B(x)','Location','northeast');xlim([0 swro_L]);ylim([swro_beta(1)*J_r/swro_x_r*.9  swro_beta(1)*J_r/swro_x_r*1.1]); hold on
-    else
-    plot(x*swro_L, swro_beta*J_r/swro_x_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('B(x)','Location','northeast');xlim([0 swro_L]); hold on
-    end
-end
-%% figure 5
-if fig(5) == 1
-figure(5); % pressure+concentration in both SWRO channels
-subplot(2,2,1);plot(x*swro_L, P_d*p_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_d(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 5;
-x2=x(2:end);C_f2=100.*C_f(2:end);
-subplot(2,2,2);plot(x*swro_L, P_f*p_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_f(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 5;
-subplot(2,2,4);plot(x2*swro_L, C_f2*C_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_f(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(2,2,3);plot(x*swro_L, 100*C_d*C_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_d(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-end
-%% figure 6
-if fig(6) == 1
-figure(6); % J_s,in and J_w,in
-x2=x(2:end);J_cross2=J_cross(2:end);J_sin2=J_sin(2:end);
-subplot(1,2,2);plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', '#77AC30', 'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{s,in}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(1,2,1);plot(x2*swro_L, J_cross2*J_r/swro_x_r, 'Color', '#0072bD', 'LineWidth',lw);xlabel('x[m]','Fontsize',10);  ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-end
-%% figure 7
-if fig(7) == 1
-f=figure(7); % PRO quantities
-f.Position = [100 1 1100 1000];
-subplot(4,3,4);plot(x*L, Q_sd*Q_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{s,d}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(4,3,1);plot(x*L, Q_wd*Q_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{w,d}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(4,3,5);plot(x*L, Q_sf*Q_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{s,f}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(4,3,2);plot(x*L, Q_wf*Q_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{w,f}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0; 
-subplot(4,3,7);plot(x*L, p_d*p_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_d^{PRO}(x)','Location','northeast'); xlim([0 L]);ay=gca; ay.YAxis.Exponent = 5;
-subplot(4,3,8);plot(x*L, p_f*p_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_f^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 5;
-subplot(4,3,10);plot(x*L, 100*c_d*C_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_d^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
-x2=x(2:end);c_f=100.*c_f(2:end);Q_cross2=Q_cross(2:end);Q_sin2=Q_sin(2:end);
-subplot(4,3,11);plot(x2*L, c_f*C_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_f^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(4,3,[3 6]);plot(x2*L, Q_cross2*Q_r/x_r, 'b', 'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(4,3,[9 12]);plot(x2*L, Q_sin2*Q_r/x_r, 'g', 'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{s,in}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
-end
-%% figure 8
-if fig(8) == 1
-f=figure(8); % SWRO quantities
-f.Position = [100 1 1100 1000];
-x2=x(2:end);J_wf2=J_wf(2:end);J_sf2=J_sf(2:end);
-subplot(4,3,1);plot(x*swro_L, J_wd*J_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{w,d}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(4,3,4);plot(x*swro_L, J_sd*J_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{s,d}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(4,3,2);plot(x2*swro_L, J_wf2*J_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{w,f}^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(4,3,5);plot(x2*swro_L, J_sf2*J_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{s,f}^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;     
-subplot(4,3,7);plot(x*swro_L, P_d*p_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_d^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 5;
-subplot(4,3,8);plot(x*swro_L, P_f*p_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_f^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 5;
-subplot(4,3,10);plot(x*swro_L, 100*C_d*C_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_d^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-C_f=100.*C_f(2:end);J_cross2=J_cross(2:end);J_sin2=J_sin(2:end);
-subplot(4,3,11);plot(x2*swro_L, C_f*C_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_f^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(4,3,[3 6]);plot(x2*swro_L, J_cross2*J_r/swro_x_r, 'Color', '#0072bD', 'LineWidth',lw);xlabel('x[m]','Fontsize',10);  ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(4,3,[9 12]);plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', '#77AC30', 'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{s,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-end
-%% figure 9
-if fig(9) == 1
 f=figure(3);  % ERD quantities
 f.Position = [100 1 900 500]; tiledlayout(1,4); 
 rot=0; rx=0; %rotation and shift of texts
@@ -628,9 +481,9 @@ lgd = legend([b(1) b(2) b(5) b(4)],'Sea water inlet','Sea water outlet','Brine o
 end
 lgd.Layout.Tile = 4; axis off ; 
 end
-%% figure 10
-if fig(9) == 1 && version(6) ==4
-close all;f=figure(3);  % ERD quantities
+%% figure 4
+if fig(3) == 1 && version(6) ==4
+close all;f=figure(4);  % ERD quantities
 f.Position = [950 1 900 500]; tiledlayout(1,4); 
 rot=0; rx=0; %rotation and shift of texts
 vals1 = [pE*p_r; pERD2*p_r; 0; p_d(end)*p_r; p_exit*p_r];
@@ -668,144 +521,12 @@ lgd = legend([b(1) b(2) b(3) b(4)],'LP ERD 2 Seawater Inlet','HP ERD 2 Seawater
 lgd.Layout.Tile = 4; axis off ; 
 end
 end
-%% figures for Paper
-if fig(11) == 1
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-f=figure(1); % SWRO quantities
-f.Position = [276.3333 196.3333 1100 1000];
-x2=x(2:end);
-% total mass flows
-subplot(3,2,1);lc='#0072bD';rc='#77AC30';
-plot(x*swro_L, J_d*J_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;xlim([0 swro_L]); hold on
-yyaxis right
-plot(x*swro_L, J_f*J_r,'Color', rc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); legend('J_{d}^{RO}(x)','J_{f}^{RO}(x)','Location','northeast');xlim([0 swro_L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-% concentrations
-subplot(3,2,2);lc='#0072bD';rc='#77AC30';
-plot(x*swro_L, 100*C_d*C_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;hold on
-yyaxis right; C_f=C_f(2:end);
-plot(x2*swro_L, 100*C_f*C_r,'Color', rc,'LineWidth',lw); ylabel('[%]','Fontsize',10); legend('C_d^{RO}(x)','C_f^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-% Pressures
-subplot(3,2,3);lc='#0072bD';rc='#77AC30';
-plot(x*swro_L, P_d*p_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 5; hold on
-yyaxis right
-plot(x*swro_L, P_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{RO}(x)','P_f^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-% Permeate flows 
-subplot(3,2,5); lc='k';rc='#b81414'; J_cross2=J_cross(2:end);J_sin2=J_sin(2:end);
-plot(x2*swro_L, J_cross2*J_r/swro_x_r, 'Color', lc, 'LineWidth',lw);xlabel('x[m]','Fontsize',10);  ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0;ylim([min( J_cross2*J_r/swro_x_r) 1.1*max( J_cross2*J_r/swro_x_r)]); hold on
-yyaxis right
-plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{RO}(x)','J_{s,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([min( J_sin2*J_r/swro_x_r) 1.1*max( J_sin2*J_r/swro_x_r)]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-% Osmotic/ Hydraulic pressure difference
-subplot(3,2,6);lc='k';rc='#b81414'; osm_diff=swro_p_osm_d(2:end)-swro_p_osm_f(2:end); a=min([min(osm_diff*p_r) min((P_d-P_f)*p_r)]); b=max( [max(osm_diff*p_r), max((P_d-P_f)*p_r)]);
-plot(x*swro_L, (P_d-P_f)*p_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ylim([a b]);ay=gca; ay.YAxis.Exponent = 6; hold on
-yyaxis right; 
-plot(x2*swro_L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{RO}(x)','\Delta \pi^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([a b]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-f=figure(2); % PRO quantities
-f.Position = [936.3333 266.3333 1.0987e+03 632];
-x2=x(2:end);
-% total mass flows
-subplot(3,2,1);lc='#0072bD';rc='#77AC30';
-plot(x*L,  Q_d*Q_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;xlim([0 swro_L]); hold on
-yyaxis right
-plot(x*L, Q_f*Q_r,'Color', rc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); legend('J_{d}^{PRO}(x)','J_{f}^{PRO}(x)','Location','northeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-% concentrations
-subplot(3,2,2);lc='#0072bD';rc='#77AC30';
-plot(x*L, 100*c_d*C_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); ay=gca; ay.YAxis(1).Exponent = 0;hold on
-yyaxis right; c_f=c_f(2:end);
-plot(x2*L, 100*c_f*C_r,'Color', rc,'LineWidth',lw); ylabel('[%]','Fontsize',10); legend('C_d^{PRO}(x)','C_f^{PRO}(x)','Location','southeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;ay.YAxis(2).Exponent = 0;
-% Pressures
-subplot(3,2,3);lc='#0072bD';rc='#77AC30';
-plot(x*L, p_d*p_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 5; hold on
-yyaxis right
-plot(x*L, p_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{PRO}(x)','P_f^{PRO}(x)','Location','northeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-% Permeate flows 
-subplot(3,2,5); lc='k';rc='#b81414'; Q_cross2=Q_cross(2:end);Q_sin2=Q_sin(2:end);
-plot(x2*L, Q_cross2*Q_r/x_r, 'Color', lc, 'LineWidth',lw);xlabel('x[m]','Fontsize',10);  ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0; hold on
-yyaxis right
-plot(x2*L, Q_sin2*Q_r/x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{PRO}(x)','J_{s,in}^{PRO}(x)','Location','northeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-% Osmotic/ Hydraulic pressure difference
-subplot(3,2,6);lc='k';rc='#b81414';osm_diff=p_osm_d(2:end)-p_osm_f(2:end); %a=min([min(osm_diff*p_r) min((p_d-p_f)*p_r)]); b=max( [max(osm_diff*p_r), max((p_d-p_f)*p_r)]);
-plot(x*L, (p_d-p_f)*p_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ay=gca;ay.YAxis.Exponent = 6; hold on; %ylim([a b]); 
-yyaxis right; 
-plot(x2*L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{PRO}(x)','\Delta \pi^{PRO}(x)','Location','northeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; %ylim([a b]);
+%ax = axes; yyaxis('left'); yyaxis('right');
+%ax.YAxis(1).Color = [0 0 0]; ax.YAxis(2).Color = [0 0 0];
 end
-%% figures for Paper - all in individual plot
-if fig(11) == 1
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-f=figure(1);% SWRO quantities
-x2=x(2:end);
-% total mass flows
-lc='#0072bD';rc='#77AC30';
-plot(x*swro_L, J_d*J_r,'Color', lc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0; xlim([0 swro_L]); hold on
-yyaxis right
-plot(x*swro_L, J_f*J_r,'Color', rc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); legend('J_{d}^{RO}(x)','J_{f}^{RO}(x)','Location','northeast');xlim([0 swro_L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_SWROa.png');
-% concentrations
-close all; f=figure(1);
-lc='#0072bD';rc='#77AC30';
-plot(x*swro_L, 100*C_d*C_r,'Color', lc,'LineWidth',lw); ylabel('[%]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;hold on
-yyaxis right; C_f=C_f(2:end);
-plot(x2*swro_L, 100*C_f*C_r,'Color', rc,'LineWidth',lw); ylabel('[%]','Fontsize',10); legend('C_d^{RO}(x)','C_f^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_SWROb.png');
-% Pressures
-close all; f=figure(6);
-lc='#0072bD';rc='#77AC30';
-plot(x*swro_L, P_d*p_r,'Color', lc,'LineWidth',lw);ylabel('[Pa]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 5; hold on
-yyaxis right
-plot(x*swro_L, P_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{RO}(x)','P_f^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_SWROc.png');
-% Permeate flows 
-close all;f=figure(1);
-lc='k';rc='#b81414'; J_cross2=J_cross(2:end);J_sin2=J_sin(2:end);
-plot(x2*swro_L, J_cross2*J_r/swro_x_r, 'Color', lc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0;ylim([min( J_cross2*J_r/swro_x_r) 1.1*max( J_cross2*J_r/swro_x_r)]); hold on
-yyaxis right
-plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{RO}(x)','J_{s,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([min( J_sin2*J_r/swro_x_r) 1.1*max( J_sin2*J_r/swro_x_r)]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_SWROd.png');
-% Osmotic/ Hydraulic pressure difference
-close all;f=figure(1);
-lc='k';rc='#b81414'; osm_diff=swro_p_osm_d(2:end)-swro_p_osm_f(2:end); a=min([min(osm_diff*p_r) min((P_d-P_f)*p_r)]); b=max( [max(osm_diff*p_r), max((P_d-P_f)*p_r)]);
-plot(x*swro_L, (P_d-P_f)*p_r,'Color', lc,'LineWidth',lw);ylabel('[Pa]','Fontsize',10); ylim([a b]);ay=gca; ay.YAxis.Exponent = 6; hold on
-yyaxis right; 
-plot(x2*swro_L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{RO}(x)','\Delta \pi^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([a b]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_SWROe.png');
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-f=figure(2); % PRO quantities
-x2=x(2:end);
-% total mass flows
-lc='#0072bD';rc='#77AC30';
-plot(x*L,  Q_d*Q_r,'Color', lc,'LineWidth',lw);ylabel('[kg/sm]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;xlim([0 swro_L]); hold on
-yyaxis right
-plot(x*L, Q_f*Q_r,'Color', rc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); legend('J_{d}^{PRO}(x)','J_{f}^{PRO}(x)','Location','northeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_PROa.png');
-% concentrations
-close all; f=figure(2);lc='#0072bD';rc='#77AC30';
-plot(x*L, 100*c_d*C_r,'Color', lc,'LineWidth',lw); ylabel('[%]','Fontsize',10); ay=gca; ay.YAxis(1).Exponent = 0;hold on
-yyaxis right; c_f=c_f(2:end);
-plot(x2*L, 100*c_f*C_r,'Color', rc,'LineWidth',lw); ylabel('[%]','Fontsize',10); legend('C_d^{PRO}(x)','C_f^{PRO}(x)','Location','southeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;ay.YAxis(2).Exponent = 0;
-set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_PROb.png');
-% Pressures
-close all; f=figure(2);lc='#0072bD';rc='#77AC30';
-plot(x*L, p_d*p_r,'Color', lc,'LineWidth',lw);ylabel('[Pa]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 5; hold on
-yyaxis right
-plot(x*L, p_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{PRO}(x)','P_f^{PRO}(x)','Location','northeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_PROc.png');
-% Permeate flows 
-close all; f=figure(2);lc='k';rc='#b81414'; Q_cross2=Q_cross(2:end);Q_sin2=Q_sin(2:end);
-plot(x2*L, Q_cross2*Q_r/x_r, 'Color', lc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0; hold on
-yyaxis right
-plot(x2*L, Q_sin2*Q_r/x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{PRO}(x)','J_{s,in}^{PRO}(x)','Location','northeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_PROd.png');
-% Osmotic/ Hydraulic pressure difference
-close all; f=figure(2);lc='k';rc='#b81414';osm_diff=p_osm_d(2:end)-p_osm_f(2:end); %a=min([min(osm_diff*p_r) min((p_d-p_f)*p_r)]); b=max( [max(osm_diff*p_r), max((p_d-p_f)*p_r)]);
-plot(x*L, (p_d-p_f)*p_r,'Color', lc,'LineWidth',lw);ylabel('[Pa]','Fontsize',10); ay=gca;ay.YAxis.Exponent = 6; hold on; %ylim([a b]); 
-yyaxis right; 
-plot(x2*L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{PRO}(x)','\Delta \pi^{PRO}(x)','Location','northeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; %ylim([a b]);
-set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_PROe.png');
 end
 
 
-%ax = axes; yyaxis('left'); yyaxis('right');
-%ax.YAxis(1).Color = [0 0 0]; ax.YAxis(2).Color = [0 0 0];
 
 
 
@@ -826,11 +547,23 @@ end
 
 
 
+%% Removed comments
 
+%Y(:,4+6)=abs(Y(:,4+6)) % ??why is this here??
 
 
+        %f_ERD=0;pMax = rho_ERD2*(ERD_eff*( p_d(end)*qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z;
+        %if pERD2 > pMax; pERD2=pMax; end
+        %if pERD2 < 0; pERD2=pMax; end  
 
+        %f_ERD=0;pMax = rho_ERD2*(ERD_eff*( p_d(1)*qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z;            
+        %if pERD2 > pMax; pERD2=pMax; end
+        %if pERD2 < 0; pERD2=pMax; end  
 
+        %f_ERD=0; pMax = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;
+        %if pERD > pMax; pERD=pMax; end
+        %if pERD < 0; pERD=pMax; end  
 
-end
-end
\ No newline at end of file
+        %f_ERD=0; pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;            
+        %if pERD > pMax; pERD=pMax; end 
+        %if pERD < 0; pERD=pMax; end  
\ No newline at end of file
diff --git a/fun_1.m b/fun_1.m
index 1cfef86..8fdab10 100644
--- a/fun_1.m
+++ b/fun_1.m
@@ -9,7 +9,7 @@ function [output1, output2, output3] = fun_1(input,option_data,obj,option_mesh,o
 %   Input:
 %       input         -   Input for the Data function
 %       option_data   -   Selects Data function used: data(input)
-%       obj           -   'SEC', 'FW', 'Rev', 'PD_net' , 'SE_net', 'approx'  or 'sol' ;
+%       obj           -   'SEC', 'FW', 'Rev', 'Pareto', 'PD_net', 'SE_net', 'SE_f', 'sol';
 %       option_mesh   -   NMax of bvp5c
 %       option_BVP    -   RelTol of bvp5c
 %
@@ -26,12 +26,17 @@ function [output1, output2, output3] = fun_1(input,option_data,obj,option_mesh,o
 %       tol=1e-6;  
 %       [Sim_01_output1, Sim_01_output2]=fun(x0,D,obj,mesh,tol),
 %
-tic;
+if nargout>2
+    tic;
+end
+
 %% Read data
 if option_data == 0; DATA = @(x)Sim_0a_data(input); end
 if option_data == 0.1; DATA = @(x)Sim_0b_data(input); end
-
-
+%
+if option_data==1; DATA = @(x)Pareto_1_data(input); end
+if option_data==2; DATA = @(x)Pareto_2_data(input); end
+%
 
 % Simulations
 if option_data == 0.11; DATA = @(x)Sim_11_data(input); end
@@ -50,13 +55,9 @@ if option_data == 0.53; DATA = @(x)Sim_PRO_03_data(input); end
 % Rough vs. Detailled
 if option_data == 0.91; DATA = @(x)Sim_rough_approx_data(input); end
 % Optimizations
-if option_data==1; DATA = @(x)Opt_1_data(input); end
-if option_data==2; DATA = @(x)Opt_2_data(input); end
 if option_data==3; DATA = @(x)Opt_3_data(input); end
 if option_data==4; DATA = @(x)Opt_4_data(input); end
 if option_data==5; DATA = @(x)Opt_L_data(input); end
-% for 3d plot
-if option_data==11; DATA = @(x)plot3d_data(input); end
 
 [H,Z,swro_Z,ro_water,ro_salt,Mw,Ms,Rw,T0,eta,sigma,p_r,rho_r,C_r,swro_L,swro_alpha,swro_R,swro_KK,swro_x_r,swro_b1,swro_b2,J_r,swro_gamma,swro_gamma2,swro_W_r,L,alpha,R,KK,x_r,b1,b2,Q_r,gamma,gamma2,W_r,cE,pE,rho_E,J_sf_0,J_wf_0,Pd_0,Pd_L,Pf_L,Q_sf_0,pd_0,pf_0,pd_L,pf_L,HP_eff,LP_eff,T_eff,V_m,ERD_eff,ERD_fric,A_ERD,eta_ERD,mix_density,pw,pe,swro_beta_fix,beta_fix,mixer_ERD,version,fig,swro_KF,swro_KD, KF, KD] ...
     = DATA(1);
@@ -114,8 +115,8 @@ if version(4) == 0 ;swro_beta=zeros(1,n); end
 % J_w,in and J_s,in
 J_cross = ((Y(:,5) - Y(:,6)) - sigma.*(swro_p_osm_d - swro_p_osm_f))./(1 + p_r*swro_alpha*swro_KK*sigma.*(swro_p_osm_d - swro_p_osm_f));
 if version(4) == 0; J_cross = ((Y(:,5) - Y(:,6)) - sigma*(swro_p_osm_d - swro_p_osm_f));  end 
-J_sin=swro_beta.*swro_del_c;
-if version(7)==1
+J_sin=swro_beta'.*swro_del_c;
+if version(7)==1 && version(4)>0
 J_cross = ((Y(:,5)-Y(:,6))-sigma.*(swro_p_osm_d - swro_p_osm_f) + swro_beta.*(Y(:,5)-Y(:,6)).*(swro_KK +1/swro_KD + 1/swro_KF)) ./ (1 + swro_alpha.*swro_beta.*(1/swro_KD + 1/swro_KF + swro_KK) + p_r*swro_alpha*sigma.*(-swro_p_osm_d./swro_KD -swro_p_osm_f.*(swro_KK + 1/swro_KF) ));
 J_sin = swro_beta.*( swro_del_c-J_cross.*( Y(:,1)./(Y(:,1) + 1)./swro_KD +  Y(:,3)./(Y(:,3) + Y(:,4)).*(swro_KK+1/swro_KF)))./(1+swro_alpha.*swro_beta.*(swro_KK+1/swro_KF+1/swro_KD));
 end
@@ -144,7 +145,7 @@ if version(5) == 0 ; beta = zeros(1,n); end
 Q_cross = (p_osm_d - p_osm_f - (Y(:,5+6) - Y(:,6+6)) - (Y(:,5+6) - Y(:,6+6)).*beta'.*KK.*alpha.*p_r)./((ones(1,n) + KK*alpha*p_r.*(beta +p_osm_f'))');
 if version(5) == 0; Q_cross = ((p_osm_d - p_osm_f) - (Y(:,5+6) - Y(:,6+6))) ; end
 Q_sin = beta'.*del_c;
-if version(8)==1 
+if version(8)==1 && version(4)>0
 Q_cross = ((p_osm_d - p_osm_f) - (Y(:,5+6) - Y(:,6+6)) - (Y(:,5+6) - Y(:,6+6)).*beta'.*(KK+1./KF+1./KD).*p_r.*alpha)./(ones(1,n) + alpha.*beta.*(KK+1./KF+1./KD) + alpha.*p_r.*(p_osm_f'.*(KK+1./KF) + p_osm_d'./KD))';
 Q_sin = beta'.*( del_c-Q_cross.*( Y(:,1+6)./(Y(:,1+6) + 1)./KD +  Y(:,3+6)./(Y(:,3+6) + Y(:,4+6)).*(KK+1/KF)))./(1+alpha.*beta');
 end
@@ -164,13 +165,14 @@ PRO_Recovery = (1- Q_f(end)./Q_f(1))*100;  % PRO recovery rate [%]
     if SWRO_Recovery < 0 
         fprintf(2,' \nERROR: Waring: SWRO recovery is negative! \n');
     end
-    if PRO_Recovery > 100 
-        fprintf(2,' \nERROR: Waring: PRO recovery is greater then 100 %% \n');
-    end
-    if PRO_Recovery < 0 
-        fprintf(2,' \nERROR: Waring: PRO recovery is negative! \n');
+    if version(6)>1
+        if PRO_Recovery > 100 
+            fprintf(2,' \nERROR: Waring: PRO recovery is greater then 100 %% \n');
+        end
+        if PRO_Recovery < 0 
+            fprintf(2,' \nERROR: Waring: PRO recovery is negative! \n');
+        end
     end
-
 %% Output of System
 J_E = (C_d(1)*J_wd(1)+J_wd(1))/2; J_ERD = J_E; J_wE = J_E/(cE+1);
 W_p1 = 1/HP_eff * (P_d(1)-pE)*(J_E *swro_Z)/rho_E; W_p1 = W_p1*swro_W_r;
@@ -192,11 +194,14 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r;
         f_2= (J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2); 
         f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
         pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_exit*Q_exit*swro_Z/rho_exit ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;
+        pMax = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_exit*Q_exit*swro_Z/rho_exit ) +  pE*J_E*swro_Z/rho_E )/J_ERD/swro_Z;
+        if pERD > pMax; pERD=pMax; end
+        if pERD < 0; pERD=pMax; end  
         W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD; W_p3=W_p3*swro_W_r; W_p4=0; W_t=0;
             if W_p3 < 0 
                 fprintf(2,' \nERROR: Waring: Pump 3 generates generates power! \n');
+                sol.stats.maxerr =1000;
             end
-
     end
     %% Version(6)=2 --> only PRO
     if version(6)==2
@@ -218,6 +223,9 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r;
                 f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2;  
                 f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
                 pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;
+                pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pE*J_E*swro_Z/rho_E)/J_ERD/swro_Z;
+                if pERD > pMax; pERD=pMax; end
+                if pERD < 0; pERD=pMax; end  
                 W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD;
                 W_t = T_eff * (p_d(end)-pE)*(Q_d(end)*Z)/local_ro_d(end); 
             end
@@ -226,6 +234,9 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r;
                 f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2  ;
                 f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
                 pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;     
+                pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pE*J_E*swro_Z/rho_E )/J_ERD/swro_Z;       
+                if pERD > pMax; pERD=pMax; end
+                if pERD < 0; pERD=pMax; end  
                 W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD;
                 W_t = T_eff * (p_d(1)-pE)*(abs(Q_d(1))*Z)/local_ro_d(1);
             end
@@ -258,13 +269,18 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r;
                 f_2= J_E*swro_Z/rho_ERD2*mix_density*(qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end)/A_ERD)^2;  
                 f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
                 pERD2 = rho_ERD2*(ERD_eff*( p_d(end)*qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z;
+                pMax = rho_ERD2*(ERD_eff*( p_d(end)*qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E )/2/J_E/swro_Z;
+                if pERD2 > pMax; pERD2=pMax; end
+                if pERD2 < 0; pERD2=pMax; end  
             end
             if version (1) == 1 %counter-current 
                 f_1= qq*Q_exit*Z/rho_exit*mix_density*(J_E_sea*swro_Z/rho_E/A_ERD)^2;  
                 f_2= J_E*swro_Z/rho_ERD2*mix_density*(qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1)/A_ERD)^2;  
                 f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
                 pERD2 = rho_ERD2*(ERD_eff*( p_d(1)*qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z;            
-
+                pMax = rho_ERD2*(ERD_eff*( p_d(1)*qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E )/2/J_E/swro_Z;            
+                if pERD2 > pMax; pERD2=pMax; end
+                if pERD2 < 0; pERD2=pMax; end 
             end     
             %calculation of Pressure P_ERD
             if version(1) == 0 %co-current
@@ -272,6 +288,9 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r;
                 f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2;  
                 f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
                 pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;
+                pMax = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pERD2*J_E*swro_Z/rho_ERD2 )/J_ERD/swro_Z;
+                if pERD > pMax; pERD=pMax; end
+                if pERD < 0; pERD=pMax; end 
                 W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD;
                 W_t = T_eff * (p_d(end)-pE)*((1-mixer_ERD)*Q_d(end)*Z)/local_ro_d(end); 
             end
@@ -279,7 +298,10 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r;
                 f_1= abs(qq*Q_d(end))*Z/local_ro_d(end)*mix_density*(J_E*swro_Z/rho_ERD2/A_ERD)^2  ;
                 f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2  ;
                 f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
-                pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;           
+                pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;
+                pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pERD2*J_E*swro_Z/rho_ERD2 )/J_ERD/swro_Z;            
+                if pERD > pMax; pERD=pMax; end 
+                if pERD < 0; pERD=pMax; end 
                 W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD;
                 W_t = T_eff * (p_d(1)-pE)*(abs((1-mixer_ERD)*Q_d(1))*Z)/local_ro_d(1);
             end
@@ -289,54 +311,52 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r;
     
 %% final output:
 W_net= - W_p1 - W_p2 - W_p3 - W_p4 + W_t; % in [W]
-PD_net = W_net./(Z*L); %in [W/m^2]    
-%old SE_net
-%if version(1)==0; SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)) + Q_d(1)*Q_r*Z/(rho_r*local_ro_d(1)))/1000/3600; end %in [kWh/m^3]
-%if version(1)==1; SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)) + abs(Q_d(end))*Q_r*Z/(rho_r*local_ro_d(end)))/1000/3600; end %in [kWh/m^3]
-%new SE_tot
-if version(1)==0; SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)))/1000/3600; end %in [kWh/m^3]
-if version(1)==1; SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)))/1000/3600; end %in [kWh/m^3]
+
+if version(6)>1
+PD_net = W_net./(Z*L); %in [W/m^2]   
+    if version(1)==0 
+        SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)) + Q_d(1)*Q_r*Z/(rho_r*local_ro_d(1)))/1000/3600; %in [kWh/m^3] 
+        SE_f   = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)))/1000/3600; %in [kWh/m^3]
+    else 
+        SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)) + abs(Q_d(end))*Q_r*Z/(rho_r*local_ro_d(end)))/1000/3600; %in [kWh/m^3] 
+        SE_f   = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)))/1000/3600; %in [kWh/m^3]
+    end 
+end
 SEC_net = W_net*(swro_local_ro_f(end)*rho_r)./(J_f(end)*J_r*swro_Z)/1000/3600; %in [kWh/m^3]  
+    if SEC_net > 0 
+        fprintf(2,' \nERROR: Waring: SEC_net is positive! \n');
+    end
 Rev= pw*FW + pe*SEC_net*FW; %in [$/h]
-tt=toc;
-
-
-
-
-
-
+if nargout>2
+    output3=toc;
+end
 
 %% if BVP-solver failed
 catch
     sol.stats.maxerr =1000;
-    SEC_net=1000;FW=0;Rev=-20;SWRO_Recovery=NaN;PRO_Recovery=NaN;PD_net=NaN;SE_net=NaN;
-    W_net=NaN;W_p1=NaN;W_p2=NaN;W_p3=NaN;W_p4=NaN;W_t=NaN;
 end
+
 %% Output of objective function
 switch (obj)
-        case 'SEC'
+        case 'SEC' % for SEC_net maximization  -->  output1 = -SEC_net
             if sol.stats.maxerr > option_BVP
-            output1 = 20;
+            output1 = 20; % using NaN causes "fmincon" to stop prematurely --> alternative "patternsearch"
             else; output1 = -SEC_net;
             end 
-        case 'FW' 
+        case 'FW' % for FW maximization  -->  output1 = -FW
             if sol.stats.maxerr > option_BVP
-            output1 = 20;
+            output1 = 20; 
             else; output1 = -FW;
             end 
-        case 'Rev' 
+        case 'Rev' % for Rev maximization  -->  output1 = -Rev
             if sol.stats.maxerr > option_BVP
             output1 = 0;
             else; output1 = -Rev;
             end
-        case 'Pareto' 
+        case 'Pareto' % SEC_net and FW maximization  -->  output1 = [-SEC_net, -FW]
             if sol.stats.maxerr > option_BVP    
-            output1 = [NaN, NaN];
+            output1 = [20, 20]; %[NaN, NaN];
             else; output1 = [-SEC_net, -FW];
-                %if SWRO_Recovery < 30;output1 = [NaN, NaN];end
-                %if SWRO_Recovery > 60;output1 = [NaN, NaN];end
-                %if PRO_Recovery < 40;output1 = [NaN, NaN];end
-                %if PRO_Recovery > 95;output1 = [NaN, NaN];end
             end 
         case 'PD_net' 
             if sol.stats.maxerr > option_BVP
@@ -348,54 +368,25 @@ switch (obj)
             output1 = 0;
             else; output1 = -SE_net;
             end
-        case 'approx' 
-            if version(6) == 0; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end
-            if version(6) == 1; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end
-            if version(6) == 2; output1=[NaN, NaN, NaN, NaN, PRO_Recovery, PD_net, SE_net]; end
-            if version(6) == 3; output1=[SEC_net, FW, Rev, SWRO_Recovery, PRO_Recovery, NaN, NaN]; end
-            if version(6) == 4; output1=[SEC_net, FW, Rev, SWRO_Recovery, PRO_Recovery, NaN, NaN]; end
-            %sol.stats.nmeshpoints
-            %sol.stats.maxerr
-            %sol.stats.nODEevals
-            %sol.stats.nBCevals
-            output2=[sol.stats.nmeshpoints sol.stats.maxerr];
-    case 'exe' 
+        case 'SE_f' 
             if sol.stats.maxerr > option_BVP
-            output1 = [NaN, NaN, NaN, NaN, NaN, NaN, NaN];
-            output2=[NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,NaN, NaN, NaN, NaN, NaN, NaN,NaN, NaN, NaN, NaN, NaN, NaN];
-            else
-            if version(6) == 0; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end
-            if version(6) == 1; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end
-            if version(6) == 2; output1=[NaN, NaN, NaN, NaN, PRO_Recovery, PD_net, SE_net]; end
-            if version(6) == 3; output1=[SEC_net, FW, SWRO_Recovery, PRO_Recovery]; end
-            if version(6) == 4; output1=[SEC_net, FW, SWRO_Recovery, PRO_Recovery]; end
-            output2=[x*swro_L; 100*C_d; J_r*J_d; p_r*P_d; p_r*P_f; x*L; 100*c_d; Q_r*Q_d; Q_r*Q_f; p_r*p_d; p_r*p_f]; % 11
-            output3=[x(2:end)*swro_L; 100*C_f(2:end); J_r*J_f(2:end); J_r/swro_x_r*J_cross(2:end); J_r/swro_x_r*J_sin(2:end); x(2:end)*L; 100*c_f(2:end); Q_r/x_r*Q_cross(2:end); Q_r/x_r*Q_sin(2:end)]; %9
-            %
-            output1=output1';output2=output2';output3=output3';
-            [a,b]=size(output1);
-            output1=[1:1:b; output1];
-            [a,b]=size(output2);
-            output2=[1:1:b; output2];
-            [a,b]=size(output3);
-            output3=[1:1:b; output3];
-            end
-    case 'sol' 
+            output1 = 0;
+            else; output1 = -SE_f;
+            end    
+        case 'sol' 
             if sol.stats.maxerr > option_BVP
-            output1 = [NaN, NaN, NaN, NaN, NaN, NaN, NaN];
-            output2=[NaN, NaN, NaN, NaN, NaN, NaN];
-            output3=NaN;
+            output1 = [NaN, NaN, NaN, NaN, NaN];
+            output2 = [NaN, NaN, NaN, NaN, NaN, NaN];
+            output3 = NaN;
             else
-            if version(6) == 0; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end
-            if version(6) == 1; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end
-            if version(6) == 2; output1=[NaN, NaN, NaN, NaN, PRO_Recovery, PD_net, SE_net]; end
-            if version(6) == 3; output1=[SEC_net, FW, Rev, SWRO_Recovery, PRO_Recovery, NaN, NaN]; end
-            if version(6) == 4; output1=[SEC_net, FW, Rev, SWRO_Recovery, PRO_Recovery, NaN, NaN]; end
+            if version(6) == 0; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN]; end
+            if version(6) == 1; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN]; end
+            if version(6) == 2; output1=[PRO_Recovery, PD_net, SE_net, SE_f, NaN]; end
+            if version(6) == 3; output1=[SEC_net, FW, Rev, SWRO_Recovery, PRO_Recovery]; end
+            if version(6) == 4; output1=[SEC_net, FW, Rev, SWRO_Recovery, PRO_Recovery]; end
             output2=[W_net, -W_p1, -W_p2,  -W_p3, -W_p4, W_t];
-            output3=tt;
             end
-  
-            
+      
 %% figure 1
 if sol.stats.maxerr ==1000
 fprintf(2,' \nERROR: bvp5c could not satisfy the relative error tolerance ---> no figures could be displayed\n');
@@ -435,13 +426,16 @@ subplot(3,2,6);lc='#0072bD';rc='#77AC30';
 plot(x*L, p_d*p_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[Pa]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 5; hold on
 yyaxis right
 plot(x*L, p_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{PRO}(x)','P_f^{PRO}(x)','Location','northeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+end
+%% figure 2
+if fig(2) == 1
 f=figure(2); %
 f.Position = [199 232.3333 1.0987e+03 639.3333];
 % Permeate flows 
 subplot(2,2,1); lc='k';rc='#b81414'; J_cross2=J_cross(2:end);J_sin2=J_sin(2:end);
 plot(x2*swro_L, J_cross2*J_r/swro_x_r, 'Color', lc, 'LineWidth',lw);xlabel('x [m]','Fontsize',10);  ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0;ylim([min( J_cross2*J_r/swro_x_r) 1.1*max( J_cross2*J_r/swro_x_r)]); hold on
 yyaxis right
-plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{RO}(x)','J_{s,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([min( J_sin2*J_r/swro_x_r) 1.1*max( J_sin2*J_r/swro_x_r)]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{RO}(x)','J_{s,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([min( J_sin2*J_r/swro_x_r) 1.1*max( J_sin2*J_r/swro_x_r)+.00001]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
 % Osmotic/ Hydraulic pressure difference
 subplot(2,2,3);lc='k';rc='#b81414'; osm_diff=swro_p_osm_d(2:end)-swro_p_osm_f(2:end); a=min([min(osm_diff*p_r) min((P_d-P_f)*p_r)]); b=max( [max(osm_diff*p_r), max((P_d-P_f)*p_r)]);
 plot(x*swro_L, (P_d-P_f)*p_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ylim([a b]);ay=gca; ay.YAxis.Exponent = 6; hold on
@@ -457,87 +451,9 @@ subplot(2,2,4);lc='k';rc='#b81414';osm_diff=p_osm_d(2:end)-p_osm_f(2:end); %a=mi
 plot(x*L, (p_d-p_f)*p_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ay=gca;ay.YAxis.Exponent = 6; hold on; %ylim([a b]); 
 yyaxis right; 
 plot(x2*L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{PRO}(x)','\Delta \pi^{PRO}(x)','Location','southeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; %ylim([a b]);
-
-
-
-end
-%% figure 2
-if fig(2) == 1
-figure(2); % pressure+concentration in both PRO channels
-subplot(2,2,1);plot(x*L, p_d*p_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('p_d(x)','Location','northeast'); xlim([0 L]);ay=gca; ay.YAxis.Exponent = 5;
-subplot(2,2,2);plot(x*L, p_f*p_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('p_f(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 5;
-subplot(2,2,3);plot(x*L, 100*c_d*C_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('c_d(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
-x2=x(2:end);c_f2=100.*c_f(2:end);
-subplot(2,2,4);plot(x2*L, c_f2*C_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('c_f(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
 end
 %% figure 3
 if fig(3) == 1
-figure(3); % Q_s,in and Q_w,in
-x2=x(2:end);Q_cross2=Q_cross(2:end);Q_sin2=Q_sin(2:end);
-subplot(1,2,1);plot(x2*L, Q_cross2*Q_r/x_r, 'b', 'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('Q_{w,in}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(1,2,2);plot(x2*L, Q_sin2*Q_r/x_r, 'g', 'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('Q_{s,in}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
-end
-%% figure 4
-lw=1.5; %Linewidth for all figures
-if fig(4) == 1
-figure(4); % Beta of PRO
-    if version(2)==0
-    plot(x*swro_L, swro_beta*J_r/swro_x_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('B(x)','Location','northeast');xlim([0 swro_L]);ylim([swro_beta(1)*J_r/swro_x_r*.9  swro_beta(1)*J_r/swro_x_r*1.1]); hold on
-    else
-    plot(x*swro_L, swro_beta*J_r/swro_x_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('B(x)','Location','northeast');xlim([0 swro_L]); hold on
-    end
-end
-%% figure 5
-if fig(5) == 1
-figure(5); % pressure+concentration in both SWRO channels
-subplot(2,2,1);plot(x*swro_L, P_d*p_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_d(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 5;
-x2=x(2:end);C_f2=100.*C_f(2:end);
-subplot(2,2,2);plot(x*swro_L, P_f*p_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_f(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 5;
-subplot(2,2,4);plot(x2*swro_L, C_f2*C_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_f(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(2,2,3);plot(x*swro_L, 100*C_d*C_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_d(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-end
-%% figure 6
-if fig(6) == 1
-figure(6); % J_s,in and J_w,in
-x2=x(2:end);J_cross2=J_cross(2:end);J_sin2=J_sin(2:end);
-subplot(1,2,2);plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', '#77AC30', 'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{s,in}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(1,2,1);plot(x2*swro_L, J_cross2*J_r/swro_x_r, 'Color', '#0072bD', 'LineWidth',lw);xlabel('x[m]','Fontsize',10);  ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-end
-%% figure 7
-if fig(7) == 1
-f=figure(7); % PRO quantities
-f.Position = [100 1 1100 1000];
-subplot(4,3,4);plot(x*L, Q_sd*Q_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{s,d}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(4,3,1);plot(x*L, Q_wd*Q_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{w,d}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(4,3,5);plot(x*L, Q_sf*Q_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{s,f}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(4,3,2);plot(x*L, Q_wf*Q_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{w,f}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0; 
-subplot(4,3,7);plot(x*L, p_d*p_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_d^{PRO}(x)','Location','northeast'); xlim([0 L]);ay=gca; ay.YAxis.Exponent = 5;
-subplot(4,3,8);plot(x*L, p_f*p_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_f^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 5;
-subplot(4,3,10);plot(x*L, 100*c_d*C_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_d^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
-x2=x(2:end);c_f=100.*c_f(2:end);Q_cross2=Q_cross(2:end);Q_sin2=Q_sin(2:end);
-subplot(4,3,11);plot(x2*L, c_f*C_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_f^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(4,3,[3 6]);plot(x2*L, Q_cross2*Q_r/x_r, 'b', 'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(4,3,[9 12]);plot(x2*L, Q_sin2*Q_r/x_r, 'g', 'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{s,in}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
-end
-%% figure 8
-if fig(8) == 1
-f=figure(8); % SWRO quantities
-f.Position = [100 1 1100 1000];
-x2=x(2:end);J_wf2=J_wf(2:end);J_sf2=J_sf(2:end);
-subplot(4,3,1);plot(x*swro_L, J_wd*J_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{w,d}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(4,3,4);plot(x*swro_L, J_sd*J_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{s,d}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(4,3,2);plot(x2*swro_L, J_wf2*J_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{w,f}^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(4,3,5);plot(x2*swro_L, J_sf2*J_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{s,f}^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;     
-subplot(4,3,7);plot(x*swro_L, P_d*p_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_d^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 5;
-subplot(4,3,8);plot(x*swro_L, P_f*p_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_f^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 5;
-subplot(4,3,10);plot(x*swro_L, 100*C_d*C_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_d^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-C_f=100.*C_f(2:end);J_cross2=J_cross(2:end);J_sin2=J_sin(2:end);
-subplot(4,3,11);plot(x2*swro_L, C_f*C_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_f^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(4,3,[3 6]);plot(x2*swro_L, J_cross2*J_r/swro_x_r, 'Color', '#0072bD', 'LineWidth',lw);xlabel('x[m]','Fontsize',10);  ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-subplot(4,3,[9 12]);plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', '#77AC30', 'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{s,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
-end
-%% figure 9
-if fig(9) == 1
 f=figure(3);  % ERD quantities
 f.Position = [100 1 900 500]; tiledlayout(1,4); 
 rot=0; rx=0; %rotation and shift of texts
@@ -577,9 +493,9 @@ lgd = legend([b(1) b(2) b(5) b(4)],'Sea water inlet','Sea water outlet','Brine o
 end
 lgd.Layout.Tile = 4; axis off ; 
 end
-%% figure 10
-if fig(9) == 1 && version(6) ==4
-close all;f=figure(3);  % ERD quantities
+%% figure 4
+if fig(3) == 1 && version(6) ==4
+close all;f=figure(4);  % ERD quantities
 f.Position = [950 1 900 500]; tiledlayout(1,4); 
 rot=0; rx=0; %rotation and shift of texts
 vals1 = [pE*p_r; pERD2*p_r; 0; p_d(end)*p_r; p_exit*p_r];
@@ -617,171 +533,8 @@ lgd = legend([b(1) b(2) b(3) b(4)],'LP ERD 2 Seawater Inlet','HP ERD 2 Seawater
 lgd.Layout.Tile = 4; axis off ; 
 end
 end
-%% figures for Paper
-if fig(11) == 1
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-f=figure(1); % SWRO quantities
-f.Position = [276.3333 196.3333 1100 1000];
-x2=x(2:end);
-% total mass flows
-subplot(3,2,1);lc='#0072bD';rc='#77AC30';
-plot(x*swro_L, J_d*J_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;xlim([0 swro_L]); hold on
-yyaxis right
-plot(x*swro_L, J_f*J_r,'Color', rc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); legend('J_{d}^{RO}(x)','J_{f}^{RO}(x)','Location','northeast');xlim([0 swro_L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-% concentrations
-subplot(3,2,2);lc='#0072bD';rc='#77AC30';
-plot(x*swro_L, 100*C_d*C_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;hold on
-yyaxis right; C_f=C_f(2:end);
-plot(x2*swro_L, 100*C_f*C_r,'Color', rc,'LineWidth',lw); ylabel('[%]','Fontsize',10); legend('C_d^{RO}(x)','C_f^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-% Pressures
-subplot(3,2,3);lc='#0072bD';rc='#77AC30';
-plot(x*swro_L, P_d*p_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 5; hold on
-yyaxis right
-plot(x*swro_L, P_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{RO}(x)','P_f^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-% Permeate flows 
-subplot(3,2,5); lc='k';rc='#b81414'; J_cross2=J_cross(2:end);J_sin2=J_sin(2:end);
-plot(x2*swro_L, J_cross2*J_r/swro_x_r, 'Color', lc, 'LineWidth',lw);xlabel('x[m]','Fontsize',10);  ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0;ylim([min( J_cross2*J_r/swro_x_r) 1.1*max( J_cross2*J_r/swro_x_r)]); hold on
-yyaxis right
-plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{RO}(x)','J_{s,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([min( J_sin2*J_r/swro_x_r) 1.1*max( J_sin2*J_r/swro_x_r)]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-% Osmotic/ Hydraulic pressure difference
-subplot(3,2,6);lc='k';rc='#b81414'; osm_diff=swro_p_osm_d(2:end)-swro_p_osm_f(2:end); a=min([min(osm_diff*p_r) min((P_d-P_f)*p_r)]); b=max( [max(osm_diff*p_r), max((P_d-P_f)*p_r)]);
-plot(x*swro_L, (P_d-P_f)*p_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ylim([a b]);ay=gca; ay.YAxis.Exponent = 6; hold on
-yyaxis right; 
-plot(x2*swro_L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{RO}(x)','\Delta \pi^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([a b]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-f=figure(2); % PRO quantities
-f.Position = [936.3333 266.3333 1.0987e+03 632];
-x2=x(2:end);
-% total mass flows
-subplot(3,2,1);lc='#0072bD';rc='#77AC30';
-plot(x*L,  Q_d*Q_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;xlim([0 swro_L]); hold on
-yyaxis right
-plot(x*L, Q_f*Q_r,'Color', rc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); legend('J_{d}^{PRO}(x)','J_{f}^{PRO}(x)','Location','northeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-% concentrations
-subplot(3,2,2);lc='#0072bD';rc='#77AC30';
-plot(x*L, 100*c_d*C_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); ay=gca; ay.YAxis(1).Exponent = 0;hold on
-yyaxis right; c_f=c_f(2:end);
-plot(x2*L, 100*c_f*C_r,'Color', rc,'LineWidth',lw); ylabel('[%]','Fontsize',10); legend('C_d^{PRO}(x)','C_f^{PRO}(x)','Location','southeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;ay.YAxis(2).Exponent = 0;
-% Pressures
-subplot(3,2,3);lc='#0072bD';rc='#77AC30';
-plot(x*L, p_d*p_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 5; hold on
-yyaxis right
-plot(x*L, p_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{PRO}(x)','P_f^{PRO}(x)','Location','northeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-% Permeate flows 
-subplot(3,2,5); lc='k';rc='#b81414'; Q_cross2=Q_cross(2:end);Q_sin2=Q_sin(2:end);
-plot(x2*L, Q_cross2*Q_r/x_r, 'Color', lc, 'LineWidth',lw);xlabel('x[m]','Fontsize',10);  ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0; hold on
-yyaxis right
-plot(x2*L, Q_sin2*Q_r/x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{PRO}(x)','J_{s,in}^{PRO}(x)','Location','northeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-% Osmotic/ Hydraulic pressure difference
-subplot(3,2,6);lc='k';rc='#b81414';osm_diff=p_osm_d(2:end)-p_osm_f(2:end); %a=min([min(osm_diff*p_r) min((p_d-p_f)*p_r)]); b=max( [max(osm_diff*p_r), max((p_d-p_f)*p_r)]);
-plot(x*L, (p_d-p_f)*p_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ay=gca;ay.YAxis.Exponent = 6; hold on; %ylim([a b]); 
-yyaxis right; 
-plot(x2*L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{PRO}(x)','\Delta \pi^{PRO}(x)','Location','northeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; %ylim([a b]);
-end
-%% figures for Paper - all in individual plot
-if fig(11) == 1
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-f=figure(1);% SWRO quantities
-x2=x(2:end);
-% total mass flows
-lc='#0072bD';rc='#77AC30';
-plot(x*swro_L, J_d*J_r,'Color', lc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0; xlim([0 swro_L]); hold on
-yyaxis right
-plot(x*swro_L, J_f*J_r,'Color', rc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); legend('J_{d}^{RO}(x)','J_{f}^{RO}(x)','Location','northeast');xlim([0 swro_L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_SWROa.png');
-% concentrations
-close all; f=figure(1);
-lc='#0072bD';rc='#77AC30';
-plot(x*swro_L, 100*C_d*C_r,'Color', lc,'LineWidth',lw); ylabel('[%]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;hold on
-yyaxis right; C_f=C_f(2:end);
-plot(x2*swro_L, 100*C_f*C_r,'Color', rc,'LineWidth',lw); ylabel('[%]','Fontsize',10); legend('C_d^{RO}(x)','C_f^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_SWROb.png');
-% Pressures
-close all; f=figure(6);
-lc='#0072bD';rc='#77AC30';
-plot(x*swro_L, P_d*p_r,'Color', lc,'LineWidth',lw);ylabel('[Pa]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 5; hold on
-yyaxis right
-plot(x*swro_L, P_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{RO}(x)','P_f^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_SWROc.png');
-% Permeate flows 
-close all;f=figure(1);
-lc='k';rc='#b81414'; J_cross2=J_cross(2:end);J_sin2=J_sin(2:end);
-plot(x2*swro_L, J_cross2*J_r/swro_x_r, 'Color', lc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0;ylim([min( J_cross2*J_r/swro_x_r) 1.1*max( J_cross2*J_r/swro_x_r)]); hold on
-yyaxis right
-plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{RO}(x)','J_{s,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([min( J_sin2*J_r/swro_x_r) 1.1*max( J_sin2*J_r/swro_x_r)]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_SWROd.png');
-% Osmotic/ Hydraulic pressure difference
-close all;f=figure(1);
-lc='k';rc='#b81414'; osm_diff=swro_p_osm_d(2:end)-swro_p_osm_f(2:end); a=min([min(osm_diff*p_r) min((P_d-P_f)*p_r)]); b=max( [max(osm_diff*p_r), max((P_d-P_f)*p_r)]);
-plot(x*swro_L, (P_d-P_f)*p_r,'Color', lc,'LineWidth',lw);ylabel('[Pa]','Fontsize',10); ylim([a b]);ay=gca; ay.YAxis.Exponent = 6; hold on
-yyaxis right; 
-plot(x2*swro_L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{RO}(x)','\Delta \pi^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([a b]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_SWROe.png');
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-f=figure(2); % PRO quantities
-x2=x(2:end);
-% total mass flows
-lc='#0072bD';rc='#77AC30';
-plot(x*L,  Q_d*Q_r,'Color', lc,'LineWidth',lw);ylabel('[kg/sm]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;xlim([0 swro_L]); hold on
-yyaxis right
-plot(x*L, Q_f*Q_r,'Color', rc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); legend('J_{d}^{PRO}(x)','J_{f}^{PRO}(x)','Location','northeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_PROa.png');
-% concentrations
-close all; f=figure(2);lc='#0072bD';rc='#77AC30';
-plot(x*L, 100*c_d*C_r,'Color', lc,'LineWidth',lw); ylabel('[%]','Fontsize',10); ay=gca; ay.YAxis(1).Exponent = 0;hold on
-yyaxis right; c_f=c_f(2:end);
-plot(x2*L, 100*c_f*C_r,'Color', rc,'LineWidth',lw); ylabel('[%]','Fontsize',10); legend('C_d^{PRO}(x)','C_f^{PRO}(x)','Location','southeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;ay.YAxis(2).Exponent = 0;
-set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_PROb.png');
-% Pressures
-close all; f=figure(2);lc='#0072bD';rc='#77AC30';
-plot(x*L, p_d*p_r,'Color', lc,'LineWidth',lw);ylabel('[Pa]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 5; hold on
-yyaxis right
-plot(x*L, p_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{PRO}(x)','P_f^{PRO}(x)','Location','northeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_PROc.png');
-% Permeate flows 
-close all; f=figure(2);lc='k';rc='#b81414'; Q_cross2=Q_cross(2:end);Q_sin2=Q_sin(2:end);
-plot(x2*L, Q_cross2*Q_r/x_r, 'Color', lc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0; hold on
-yyaxis right
-plot(x2*L, Q_sin2*Q_r/x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{PRO}(x)','J_{s,in}^{PRO}(x)','Location','northeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
-set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_PROd.png');
-% Osmotic/ Hydraulic pressure difference
-close all; f=figure(2);lc='k';rc='#b81414';osm_diff=p_osm_d(2:end)-p_osm_f(2:end); %a=min([min(osm_diff*p_r) min((p_d-p_f)*p_r)]); b=max( [max(osm_diff*p_r), max((p_d-p_f)*p_r)]);
-plot(x*L, (p_d-p_f)*p_r,'Color', lc,'LineWidth',lw);ylabel('[Pa]','Fontsize',10); ay=gca;ay.YAxis.Exponent = 6; hold on; %ylim([a b]); 
-yyaxis right; 
-plot(x2*L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{PRO}(x)','\Delta \pi^{PRO}(x)','Location','northeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; %ylim([a b]);
-set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_PROe.png');
-end
-
-
 %ax = axes; yyaxis('left'); yyaxis('right');
 %ax.YAxis(1).Color = [0 0 0]; ax.YAxis(2).Color = [0 0 0];
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
 end
 end
 
@@ -808,50 +561,4 @@ end
 
 %% Removed comments
 
-%Y(:,4+6)=abs(Y(:,4+6)) % ??why is this here??
-
-% -3.9191    0.3247    0.0131   32.8727    7.1201       NaN       NaN
-%
-% 1.0e+03 *
-% -1.2742   -0.9744         0   -0.4873   -0.0010    0.1885
-
-% In ideal cases take mean of some quantities. 
-% In ideal cases salt will fluctuate around 1e-16. We take mean here for 
-% 
-% This is just for figures. (), otherwise one could use ylim ) 
-%if version(4)==0 ; J_f=J_wf; J_sf=zeros(1,n); J_sd=mean(J_sd)*ones(1,n); C_f=zeros(1,n); end
-%if version(5)==0 ; Q_f=Q_wf; Q_sf=zeros(1,n); Q_sd=mean(Q_sd)*ones(1,n); c_f=zeros(1,n); end
-
-        %f_ERD=0; pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_exit*Q_exit*swro_Z/rho_exit ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;
-        %if pERD > pMax; pERD=pMax; end
-        %if pERD < 0; pERD=pMax; end  
-        %display(['pERD:              ',num2str(pERD)]) % pERD use in the model
-        %display(['pERD_theoretical:  ',num2str(pMax)]) % theoretical pERD with no leak, friction perm and 100% efficiency
-
-        %f_ERD=0;pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;
-        %if pERD > pMax; pERD=pMax; end
-        %if pERD < 0; pERD=pMax; end  
-
-        %f_ERD=0;pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;       
-        %display(['pERD:              ',num2str(pERD)]) % pERD use in the model
-        %display(['pERD_theoretical:  ',num2str(pMax)]) % theoretical pERD with no leak, friction perm and 100% efficiency
-        %if pERD > pMax; pERD=pMax; end
-        %if pERD < 0; pERD=pMax; end     
-
-        %J_wERD = J_E/(C_ERD+1);
-
-        %f_ERD=0;pMax = rho_ERD2*(ERD_eff*( p_d(end)*qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z;
-        %if pERD2 > pMax; pERD2=pMax; end
-        %if pERD2 < 0; pERD2=pMax; end  
-
-        %f_ERD=0;pMax = rho_ERD2*(ERD_eff*( p_d(1)*qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z;            
-        %if pERD2 > pMax; pERD2=pMax; end
-        %if pERD2 < 0; pERD2=pMax; end  
-
-        %f_ERD=0; pMax = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;
-        %if pERD > pMax; pERD=pMax; end
-        %if pERD < 0; pERD=pMax; end  
-
-        %f_ERD=0; pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;            
-        %if pERD > pMax; pERD=pMax; end 
-        %if pERD < 0; pERD=pMax; end  
\ No newline at end of file
+%Y(:,4+6)=abs(Y(:,4+6)) % ??why is this here?? 
\ No newline at end of file
diff --git a/nonlcon.m b/nonlcon.m
index e9177f0..b21be8c 100644
--- a/nonlcon.m
+++ b/nonlcon.m
@@ -1,4 +1,4 @@
-function [c, ceq] = nonlcon(input,option,option_data,option_mesh,option_BVP, minFW)
+function [c, ceq] = nonlcon(input,option,option_data,option_mesh,option_BVP, minConstraint)
 %%  nonlcon  Nonlinear condition in fmincon 
 %
 %       nonlcon(input,option,option_data,option_mesh,option_BVP, minFW)
@@ -13,7 +13,7 @@ function [c, ceq] = nonlcon(input,option,option_data,option_mesh,option_BVP, min
 %       option_data   -   Selects Data function used: data(input)
 %       option_mesh   -   NMax of bvp5c
 %       option_BVP    -   RelTol of bvp5c
-%       minFW         -   lower bound for FW
+%       minConstraint -   lower bound for FW/SEC_net
 %
 %   Output:
 %       c(x) ≤0
@@ -31,7 +31,10 @@ switch (option)
     c=[];
     case 'FW'
     % output from fun: -FW    
-    c= fun_1(input,option_data,'FW',option_mesh,option_BVP) + minFW; 
+    c= fun_1(input,option_data,'FW',option_mesh,option_BVP) + minConstraint; 
+    case 'SEC'
+    % output from fun: -SEC_net    
+    c= fun_1(input,option_data,'SEC',option_mesh,option_BVP) + minConstraint; 
 end
 ceq=[];
 end
\ No newline at end of file
-- 
GitLab