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

clc
clear all
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.0; %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; 
v_3_nom = - A * W_b4_nom;

x0 = [h_2; W_a4_nom; W_b4_nom];

%% 2. Disturbances

q_1tstep = 1000; 
q_2tstep = 3000;

q_1step = -2;
q_2step = 0.65;

%% 3. Controller Parameters


k_2 = 1;
tauc_2 = 500;
theta_2 = 0;
tau_2 = 24.2;

Kc_2 = (1/k_2)*(tau_2/(tauc_2+theta_2));
tauI_2 = min(tau_2,4*(tauc_2+theta));


k_3 = -1/A;                                         %gain of y(s)/v(s)[min]
tauc_3 = 105;                                              %tuning parameter
theta_3 = 0;                                                    %time delay
tau_3 = -1/A;

Kc_3 = (1/k_3)*(tau_3/(tauc_3+theta_3));
tauI_3 = min(tau_3,4*(tauc_3+theta_3));

k_4 = -796.77;
tauc_4 = 0.0005;

KI_4 = 1/(k_4*tauc_4);

k_5 = 0.8177;
tauc_5 = 0.0005;

KI_5 = 1/(k_5*tauc_5);
%% 4. Run

sim('CL_model_SecondB.slx');

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

subplot(121)
plot(t,pH)
set(gca,'TickLabelInterpreter','latex')
xlabel('$t$ [s]','Interpreter','latex')
ylabel('pH','Interpreter','latex')

subplot(122)
plot(t,q_3_GLT)
set(gca,'TickLabelInterpreter','latex')
xlabel('$t$ [s]','Interpreter','latex')
ylabel('$q_3$ [ml/s]','Interpreter','latex')

figure()
subplot(121)
plot(t,Wb4sp,t,W_b4)
ylabel('W_{b4}')
xlabel('t')

subplot(122)
plot(t,W_a4)


% filename = 'Alt_Wb4';
% save(filename,'y_sp','y','x_sp','Wb4sp','v2sp','t','x','pH',...
%        'H_plus','q_3_GLT','W_a4','W_b4','in1','in2','v3')
