%% This script is used to predict pH transients
% Equations adapted from Pabst and Carta 2006 using Eqs 34-36. This script
% is designed to be used for a monovalent buffer where the parameter values
% are for sodium acetate.

clear; close all; tic;


%% Inputs

L = 2.5;                         % Column length (cm)
D = 0.5;                         % Column diameter (cm)
F = 1;                           % Flow rate (CV/min)
Ee = 0.4;                        % Interstitial porosity
Ep = 0.91;                       % Intraparticle porosity
Et = Ee + (1-Ee)*Ep;             % Total porosity
pKa = 4.76;                      % pKa of buffer species
pK = 4.64;                       % pKa of ligand
qR = 0.08;                       % Ligand density (M packed column)
B = 0.0172;                      % Na+ concentration in buffer (M)
S_i = 0;                         % Initial NaCl concentration (M)
S = 0.5;                         % NaCl concentration in elution buffer (M)
A = 0.5114;                      % Temperature effects constant (298 K)
b = 0.16;                        % Buffer parameter
zA = 0;                          % Net charge of conjugate acid species
zI = 1;                          % Net charge of H+
pHi = 5.5;                       % Initial pH
N = 60;                          % Number of axial discretizations
k_kin = 1e8;                     % Kinetic rate for LDF

% Salt sequence (start and end salt for each section), four sections
% intended for load, wash, elution, strip
csalt0 = [B+S_i B+S_i B+S_i+S B+S_i+S];
csalt1 = [B+S_i B+S_i B+S_i+S B+S_i+S];
cv_seq = [2 3 10 5];


%% Preliminary calculations

Ac = 0.25*pi*D^2;                % Cross sectional area (cm^2)
V = Ac*L;                        % Column volume (ml)
phi = (1-Ee)/Ee;                 % Phase ratio
I = B + S;                       % Ionic strength (mM)
u = F*V/Ac;                      % Linear velocity (cm/min)
nu = u/Ee;                       % Interstitial velocity (cm/min)
K = 10^-pK;                      % Ka of ligand
qR = qR/(1-Et);                  % Convert to resin volume


%% Region I calculations

cNaI = B + S_i;                  % Sodium concentration at region I
cClI = S_i;                      % Chloride concentration at region I
gI = ...                         % Activity coefficient of H+ at region I
    10^(-zI^2*(A*sqrt(cNaI)/(1+sqrt(cNaI)) - b*cNaI)); 
pKaI = ...                       % pKa of buffer at region III
    pKa + 2*(zA-1)*(A*sqrt(cNaI)/(1+sqrt(cNaI)) - b*cNaI);
cHI = 10^-pHi/gI;                % H+ concentration at region I
cAI = ...                        % Buffer concentration at region I
    (cHI/(10^-pKaI) + 1)*(cNaI-cClI);   
qNaI = ...                       % Solid phase conc. at region I
    0.5*(-K*cNaI/cHI + sqrt((K*cNaI/cHI)^2 + 4*qR*K*cNaI/cHI));


%% Prepare arrays for ODE solver 

% Calculate time series for gradient steps                             
t_seq = cv_seq/F*60;
t_factor = 1;
t_sim = linspace(0,sum(t_seq),(ceil(sum(t_seq))+1)*t_factor);

% Initialize arrays and apply initial conditions
cA = ones(N+1,1)*cAI;              % Acetate
cCl = ones(N+1,1)*cClI;            % Chloride
cNa = ones(N+1,1)*cNaI;            % Liquid phase Na
qNa = ones(N+1,1)*qNaI;            % Solid phase Na
qNa(1) = 0;
y = [cA; cCl; cNa; qNa];

% Initialize model structure for holding inputs
ms.pKa = pKa; ms.zA = zA; ms.A = A; ms.b = b; ms.K = K; ms.qR = qR; 
ms.S = S; ms.L = L; ms.D = D; ms.u = u; ms.Et = Et;
ms.csalt0 = csalt0; ms.csalt1 = csalt1; ms.t_seq = t_seq; ms.N = N;
ms.pHi = pHi; ms.k_kin = k_kin; ms.zI = zI; ms.zA = zA; ms.B = B;

% Solving the model
options = odeset('RelTol',1e-6,'AbsTol',1e-8,'NormControl','on');
% options = odeset('NormControl','on');
[t, y]  =  ode15s(@dsys,[0 sum(t_seq)],y,options,ms);   
% odeset('InitialStep',1,'MaxStep',1)

% Parse column profiles
cA = y(:,0*N+1:1*N+1);          % Acetate
cCl = y(:,1*N+2:2*N+2);         % Chloride
cNa = y(:,2*N+3:3*N+3);         % Liquid phase Na
qNa = y(:,3*N+4:4*N+4);         % Solid phase Na

% Calculate parameters and pH
Ka = 10.^-(pKa + 2*(zA-1)*(A*sqrt(cNa)./(1+sqrt(cNa))-b*cNa));
g = 10.^-(zI^2*(A*sqrt(cNa)./(1+sqrt(cNa))-b*cNa));
cH = Ka.*(cA./(cNa-cCl) - 1);
pH = -log10(g.*cH);

% Assign outlet profiles
cA_out = cA(:,end);
cCl_out = cCl(:,end);
cNa_out = cNa(:,end);
qNa = qNa(:,end);
pH_out = pH(:,end);
CV_out = t*(F/60);

% Generate plot with pH and salt
figure(1);
yyaxis left
plot(CV_out,pH_out,'k','linewidth',2)
set(gca,'xminortick','on','yminortick','on','fontsize',12,'ycolor','k')
ylim([4 7])
ylabel('pH')
yyaxis right
plot(CV_out,cNa_out,'color',[0.93 0 0],'linewidth',1.5)
ylabel('[Na+] (M)')
ylim([0 0.6])
xlabel('CV')
set(gca,'xminortick','on','yminortick','on','fontsize',12,'ycolor',[0.93 0 0])
grid on
box on

toc;

function dydt = dsys(t,y,ms)
    % Parse inputs from model structure
    t_seq = [0 ms.t_seq]; csalt0 = ms.csalt0; csalt1 = ms.csalt1;
    pHi = ms.pHi; A = ms.A; b = ms.b; K = ms.K; qR = ms.qR;
    Et = ms.Et; u = ms.u/60; k_kin = ms.k_kin; N = ms.N; h = ms.L/N;
    zI = ms.zI; zA = ms.zA; pKa = ms.pKa; B = ms.B;
    
    % Reading from y vector to single vectors
    cA(:,1) = y(0*N+1:1*N+1);
    cCl(:,1) = y(1*N+2:2*N+2);
    cNa(:,1) = y(2*N+3:3*N+3);
    qNa(:,1) = y(3*N+4:4*N+4);
    
    % Zeros dot vectors
    dydt = zeros(4*N+4,1); dcAdx = zeros(N+1,1); dcCldt = zeros(N+1,1);
    dcNadt = zeros(N+1,1); dqNadt = zeros(N+1,1); cH = zeros(N+1,1);
    pH = zeros(N+1,1); g = zeros(N+1,1); Ka = zeros(N+1,1);

    % Boundary conditions at inlet (for salt gradient)
    for i=1:numel(t_seq)-1
        if t > sum(t_seq(1:i)) && t <= sum(t_seq(1:i+1))
            cNa(1) = csalt0(i) + ((csalt1(i)-csalt0(i))/t_seq(i+1))*(t-sum(t_seq(1:i)));
            cCl(1) = csalt0(i) + ((csalt1(i)-csalt0(i))/t_seq(i+1))*(t-sum(t_seq(1:i))) - B;            
        end
    end
    
    % Boundary conditions at inlet
    pH(1) = pHi; 
    g(1) = 10^-(zI^2*(A*sqrt(cNa(1))/(1+sqrt(cNa(1)))-b*cNa(1)));
    cH(1) = 10^-pH(1)/g(1); 
    Ka(1) = 10^-(pKa + 2*(zA-1)*(A*sqrt(cNa(1))/(1+sqrt(cNa(1)))-b*cNa(1)));
    cA(1) = (cH(1)/Ka(1) + 1)*(cNa(1) - cCl(1));
    
    % Set up differential balances for all discretization cells
    for i=2:N+1
        dcAdx(i) = u/Et*(cA(i-1) - cA(i))/h;
        dcCldt(i) = u/Et*(cCl(i-1) - cCl(i))/h;
        Ka(i) = 10^-(pKa + 2*(zA-1)*(A*sqrt(cNa(i))/(1+sqrt(cNa(i)))-b*cNa(i)));
        g(i) = 10^-(zI^2*(A*sqrt(cNa(i))/(1+sqrt(cNa(i)))-b*cNa(i)));
        cH(i) = Ka(i)*(cA(i)/(cNa(i)-cCl(i)) - 1);
        pH(i) = -log10(g(i)*cH(i));
%         qNa_eq = 0.5*(-K*cNa(i)/cH(i) + ...
%             sqrt((K*cNa(i)/cH(i))^2-4*qR*K*cNa(i)/cH(i)));
%         dqNadt(i) = k_kin*(qNa(i) - qNa_eq);
        cNa_eq = qNa(i)^2*cH(i)/((qR-qNa(i))*K);
        dqNadt(i) = k_kin*(cNa(i) - cNa_eq);
        dcNadt(i) = -(1-Et)/Et*dqNadt(i) + u/Et*(cNa(i-1) - cNa(i))/h;
    end

    % Transfer back to ydot vector for solving
    dydt(0*N+1:1*N+1) = dcAdx; dydt(1*N+2:2*N+2) = dcCldt;
    dydt(2*N+3:3*N+3) = dcNadt; dydt(3*N+4:4*N+4) = dqNadt;
end



