clc
clear all
tic
%F = 0.001:1:1600;
F = 50:10:1600;
V = 10;
Tk = 128.95 + 273.15;
y0 = [1 4 Tk Tk];
Cain = 5.1;
Qk = -4495.7;
Kw = 4032;
Ar = 0.215;

rpar(1)=Cain;
rpar(2)=Qk;

op=optimset('Display','off');

for k = 1:length(F)
    rpar(3) = F(k);
    rpar(1) = 4;
    estados_est(k,:) = fsolve('vanvusse_EE',y0,op,rpar);

    rpar(1) = 5.1;
    estados_est2(k,:) = fsolve('vanvusse_EE',y0,op,rpar);
    
    rpar(1) = 6;    
    estados_est3(k,:) = fsolve('vanvusse_EE',y0,op,rpar);
end
F_V = F/V;

figure(2)
plot(F_V,estados_est(:,2),'-.',F_V,estados_est2(:,2),'-',F_V,estados_est3(:,2),':')
xlabel('F/V (h^-^1)')
ylabel('Cb_e_e (mol/L)')
legend('Ca_{in}=4.0 mol/L','Ca_{in}=5.1 mol/L','Ca_{in}=6.0 mol/L');

figure(3)
plot(F_V,estados_est(:,3),'-.',F_V,estados_est2(:,3),'-',F_V,estados_est3(:,3),':')
xlabel('F/V (h^-^1)')
ylabel('T_e_e (mol/L)')
legend('Ca_{in}=4.0 mol/L','Ca_{in}=5.1 mol/L','Ca_{in}=6.0 mol/L');

rpar(1)=Cain;
rpar(3)=700;
Q = -5000*[4:-0.1:-4];

for k = 1:length(Q)
    rpar(2) = Q(k);
    rpar(1) = 4;
    estados_est4(k,:) = fsolve('vanvusse_EE',y0,op,rpar);

    rpar(1) = 5.1;
    estados_est5(k,:) = fsolve('vanvusse_EE',y0,op,rpar);
    
    rpar(1) = 6;    
    estados_est6(k,:) = fsolve('vanvusse_EE',y0,op,rpar);
end
QK=-Q/(Kw*Ar);

figure(4)
plot(QK,estados_est4(:,2),'-.',QK,estados_est5(:,2),'-',QK,estados_est6(:,2),':')
xlabel('-Q/(K_w A_r) (^oC/h)')
ylabel('Cb_e_e (mol/L)')
legend('Ca_{in}=4.0 mol/L','Ca_{in}=5.1 mol/L','Ca_{in}=6.0 mol/L');

figure(5)
plot(QK,estados_est4(:,3),'-.',QK,estados_est5(:,3),'-',QK,estados_est6(:,3),':')
xlabel('-Q/(K_w A_r) (^oC/h)')
ylabel('T_e_e (mol/L)')
legend('Ca_{in}=4.0 mol/L','Ca_{in}=5.1 mol/L','Ca_{in}=6.0 mol/L');

toc
