%pH Neutralisation
%PI Only
%CK
%20/10/2020
% u = q3; y = pH in q4/tank


close all
set(0,'DefaultTextFontName','Times',...
     'DefaultTextFontSize',20,...
     'DefaultAxesFontName','Times',...
     'DefaultAxesFontSize',20,...
     'DefaultLineLineWidth',1,...
     'DefaultLineMarkerSize',7.75,...
     'DefaultStairLineWidth',2)

%% 1. Initialisation

%parameters

q1_conc = 0.003; %1[M HNO3]
q2_conc = 0.03; %2[NaHCO3]
q3_conc_NaOH = 0.003; %3[M NaOH]
q3_conc_NaHCO3 = 0.0005; %4[NaHCO3]
A_1 = 207; %5 [cm^2]
A_2 = 43; %6[cm^2]
z = 11.5; %7[cm]
n = 0.607;%8 [-]
K_a1 = 4.47e-7; %9[-]
K_a2 = 5.62e-11; %10[-]
K_w = 1.00e-14; %11[-]
del_t = 15.0; %12[s]
del_tc = 1.0;%13[s]
tau_pH = 15.0; %14[s]
tau_h = 15.0; %15[s]
tau_v = 6.0; %16[s]
theta = 10.0; %17[s]
q_1 = 16.6; %18[ml/s]
q_2 = 0.55; %19[ml/s]
q_3 = 15.550262138347007; %20[ml/s]
q_1e = 16.6; %21[ml/s]
q_4 = 32.8; %22[ml/s]
W_a1 = 3e-3; %23[M]
W_b1 = 0; %24[M]
W_a2 = -0.03; %25[M}
W_b2 = 0.03; %26[M]
W_a3 = -3.05e-3; %27[M]
W_b3 = 5e-5; %28[M]
h_1 = 14.0; %29[cm]
h_2 = 3.0; %30[cm]
W_a4_nom =-4.320546258065119e-04;% -4.320336391437310e-04; %31[M]
W_b4_nom =5.283600796171088e-04; %5.283639143730888e-04; %32[M]
pH_nom = 7.000000; %7.026296087741869; %33[-]
%A =(W_b3-W_b4_nom)/(A_1*h_1); %34
A = (-q_1-q_2-q_3)/(A_1*h_1);


p = [q1_conc;q2_conc;q3_conc_NaOH;q3_conc_NaHCO3;A_1;A_2;z;n;K_a1;K_a2;...
    K_w;del_t;del_tc;tau_pH;tau_h;tau_v;theta;q_1;q_2;q_3;q_1e;q_4;...
    W_a1;W_b1;W_a2;W_b2;W_a3;W_b3;h_1;h_2;W_a4_nom;W_b4_nom;pH_nom; A];

H_plus_nom = 10^(-pH_nom);
v_2_nom = H_plus_nom; 
c1_nom = (K_a1/H_plus_nom + (2*K_a1*K_a2)/H_plus_nom^2)/...
    (1+K_a1/H_plus_nom+K_a1*K_a2/H_plus_nom^2);
W_4_nom = W_a4_nom+c1_nom*W_b4_nom;
v_3_nom = W_4_nom;

x0 = [h_2; W_a4_nom; W_b4_nom];

%% 2. Disturbances

q_1tstep = 500; 
q_2tstep = 500;

q_1step = 0;
q_2step = 0;

pH_newsp =  6.5; 

v_3_newsp = 10^(-pH_newsp)-K_w/(10^(-pH_newsp));

%% 3. Controller Parameters

k_3 = 0.002122;                                                       %gain
tauc_3 = 75;                                              %tuning parameter
theta_3 = 0;                                                    %time delay
tau_3= 150;                                                  %time constant

Kc_3 = (1/k_3)*(tau_3/(tauc_3+theta_3));         %Inner Controller Kc(SIMC)
tauI_3 = min(tau_3,4*(tauc_3+theta_3));      %Inner Controller tau_I (SIMC)

k_4 = -796.77;
tauc_4 = 0.05;

KI_4 = 1/(k_4*tauc_4);

%% 4. Run

inner_controller = 1;                   % =0 for inner off, =1 for inner on
c_1_choice = 1;             % =0 to use [H+] in c1, =1 to use [H+]sp for c1

sim('model_PI.slx');

h_2 = states(:,1);
W_a4 = states(:,2);
W_b4 = states(:,3);

% c_1_choice = 0;
% sim('CL_model_combination_alt.slx');
% 
% h_2_v = states(:,1);
% W_a4_v = states(:,2);
% W_b4_v = states(:,3);
% t_v =t;
% pH_v = pH;
% q_3_v = q_3_GLT;
% v_v = v;
% H_plus_v = H_plus;
% 
% 
% inner_controller = 1; 
% outer_controller = 0;
% sim('CL_model_combination_alt.slx')
% h_2_i = states(:,1);
% W_a4_i = states(:,2);
% W_b4_i = states(:,3);
% t_i =t;
% pH_i = pH;
% q_3_i = q_3_GLT;
% v_i = v;
% H_plus_i = H_plus;
% 
% inner_controller = 1;
% outer_controller = 1;
% sim('CL_model_combination_alt.slx')
% h_2_o = states(:,1);
% W_a4_o = states(:,2);
% W_b4_o = states(:,3);
% t_o =t;
% pH_o = pH;
% q_3_o = q_3_GLT;
% v_o = v;
% H_plus_o = H_plus;
% 
% subplot(121)
% plot(t_v_sp,pH_v_sp,t_v,pH_v,t_i,pH_i,t_o,pH_o)
% set(gca,'TickLabelInterpreter','latex')
% xlabel('$t$ [s]','Interpreter','latex')
% ylabel('pH','Interpreter','latex')
% legend('No Controller, $c=f([H^+_{sp}])$','No Controller, $c=f([H^+])$','Inner Controller Only',...
%     'Outer + Inner Controller','Interpreter','latex')
% 
% subplot(122)
% plot(t_v_sp,q_3_v_sp,t_v,q_3_v,t_i,q_3_i,t_o,q_3_o)
% set(gca,'TickLabelInterpreter','latex')
% xlabel('$t$ [s]','Interpreter','latex')
% ylabel('$q_3$ [ml/s]','Interpreter','latex')
% 
% figure()
% plot(t_v_sp,v_v_sp,t_v,v_v,'--',t_i,v_i,t_o,v_o)
% set(gca,'TickLabelInterpreter','latex')
% xlabel('$t$ [s]','Interpreter','latex')
% ylabel('$v$ [M]','Interpreter','latex')
% legend('No Controller, $c=f([H^+_{sp}])$','No Controller, $c=f([H^+])$','Inner Controller Only',...
%     'Outer + Inner Controller','Interpreter','latex')
% 
% 
% % figure()
% % plot(H_plus,func_H_plus)
% % xlabel('[H^+]')
% % ylabel('f([H^+])')
% 
% 
%  filename = 'trial_pos10%';
% save(filename,'ysp','t','x','pH',...
%        'H_plus','q_3_GLT','W_a4','W_b4')
% % % 
% 
% 
plot(t,pH)
yyaxis right
plot(t,H_plus)
