% Simulates dynamic CSTR for A-to-B reaction
clc
clear

% Specify parameters
    % Scaling parameters
    gamma = -2.0; % Linear scaling slope, (-)
    delta = -1.4; % Equivalent surface enthalpy, (eV) ** enter as negative of value listed in manuscript
    alpha = 0.6; % BEP slope, (-)
    beta = 102; % BEP offset, (kJ/mol)

% Waveform parameters
    fosc = 1e0; % Frequency, (Hz)
    delU = 0.4; % Amplitude, (eV)
    BEa0 = 1.41; % Waveform centerpoint, (eV)



% Generate simulation results
    qsdot = 0.05; % mL/min, initial flowrate guess
    tic
    [SV,tos,TOFbavg,Theta_avg,qsve,PlotData] = dynamic_catalysis_CSTR(gamma,delta,alpha,beta,fosc,delU,BEa0,qsdot);
    results = [delU,fosc,TOFbavg,Theta_avg,SV,tos]
    toc
    
    
    % Variable names
    %     TOFbavg = time-averaged, dynamic turnover frequency of B (1/s)
    %     Theata_avg [Theta_A*, Theta_B*, Theta_*] = time-averaged, dynamic steady-state surface coverages of A, B, empty sites (-)
    %     PlotData = [tsx,xsv,BEasv]
    %         tsv = time on stream (s)
    %         xsv = [A(g), B(g), ThetaA*, ThetaB*] on-stream data (M or -)
    %         BEasv = BEa on-stream data (eV)
    %         TOFbinst = TOFb on-stream data (1/s)
    %     SV = space velocity at 1% time-averaged conversion of B (1/s)
    %     tos = time on stream (s)



%% Functions
function [SV,tos,TOFbavg,Theta_avg,qsve,PlotData] = dynamic_catalysis_CSTR(gamma,delta,alpha,beta,fosc,delU,BEa0,qsdot)
% Specify parameters
    Tc = 200; % reactor temperature (deg C)

% Constants:
    % Gas constant
    R = 8.31446261815324; % J/mol-K
    Rg = R*0.01; % L-bar/mol-K

    % Conditions:
    % Reactor temperature
    T = Tc + 273.15; % K

    % Feed pressure (bar)
    Pf(1) = 100; % A(g)
    Pf(2) = 0; % B(g)

    % Feed concentration
    Cf = Pf/Rg/T; % M

% Initial conditions
    x0 = [Cf';0;0];
        
% Volumetric flowrate
qdot = qsdot/60000; % L/s

% Catalyst weight
w = 138; % mg
% Catalyst density
rho = 0.53; % g/mL
% Reactor bed volume (L)
V = w*0.001/rho*0.001;

% Space velocity
SV = qdot/V; % 1/s

% Static parameters:
% Vapor phase enthalpy (kJ/mol)
H(1) = 0; % A(g)
H(2) = 0; % B(g)


% Set endpoints based on delU and centedrpoint
UR = delU/2;
UL = UR - delU;


% Oscillation period
tosc = 1/fosc; % s
tau(1) = 1/fosc/2; % Strong binding
tau(2) = tosc - tau(1); % Weak binding

% Number of oscillations
Nosc = max(fosc,20);

% Compute rate constants (1/s):
k(1,:) = cstr_constants(BEa0,UR,gamma,H,delta,alpha,beta,T);
k(2,:) = cstr_constants(BEa0,UL,gamma,H,delta,alpha,beta,T);

% Target conversion
C = 1.0; % percent

% ss tolerance
ssTol = 0.0001;

% Calculations:
    % Simulate square waveform
    [SV,tos,TOFbavg,Theta_avg,qsve,PlotData] = square_waveform_CSTR(Nosc,k,T,SV,Cf,tau,x0,C,qdot,UL,UR,BEa0,ssTol);
end

function k = cstr_constants(BEa0,delBEa,gamma,H,delta,alpha,beta,T)
% Constants:
% Gas constant
R = 8.31446261815324; % J/mol-K
R = R*0.001; % kJ/mol-K

% Boltzmann constant
kB = 1.380649e-23; % J/K

% Planck constant
h = 6.62607015e-34; % J-s

% Parameters:
% Binding energy (eV)
BE(1) = BEa0 + delBEa; % A*
BE(2) = gamma*BE(1) + H(2)/96.485 - gamma*H(1)/96.485 - (1 - gamma)*delta; % B*
% Restrict to positive values
BE = max(BE,0);

% Heat of reaction for each elementary step (kJ/mol)
delH(1) = -BE(1)*96.485; % A(g) + * <--> A*
delH(2) = (H(2) - BE(2)*96.485) - (H(1) - BE(1)*96.485); % A* <--> B*
delH(3) = BE(2)*96.485; % B* <--> B(g) + *

% Activation energy (kJ/mol)
Ea = alpha*delH(2) + beta; % BEP
% Restrict to positive values
Ea = max(Ea,delH(2));
Ea = max(Ea,0);

% Entropy of reaction (J/mol-K)
delS(1) = -135; % A(g) + * <--> A*
delS(2) = 0; % A* <--> B*
delS(3) = 135; % B* <--> B(g) + *

% Gibbs free energy of reaction 
delG = delH - T*delS*0.001; % 

% Equilibrium constants
K = exp(-delG/R/T); % unitless

% Compute rate constants:
% Forward (1/s)
k(1) = kB*T/h*exp(delS(1)*0.001/R); % A(g) + * --> A*
k(3) = kB*T/h*exp(-Ea/R/T); % A* --> B*
k(6) = kB*T/h*exp(-delS(3)*0.001/R); % B(g) + * --> B*

% Reverse (1/s)
k(2:2:4) = k(1:2:3)./K(1:2); 
k(5) = k(6)*K(3); % B* --> B(g) + *
end

function xdot = cstr(t,x,k,T,SV,Cf)
% Constants:
% Gas constant
R = 8.31446261815324; % J/mol-K
Rg = R*0.01; % L-bar/mol-K

% Reference pressure
Po = 1.01325; % bar

% Parameters:
% Catalyst weight
w = 138; % mg
% Active site density
srho = 20; % micromol/g
% Number of active sites (gmol)
Nsites = w*0.001*srho*0.000001;

% Catalyst density
rho = 0.53; % g/mL
% Reactor bed volume (L)
V = w*0.001/rho*0.001;

% Stoichiometry matrix (unitless)
nu(1,:) = [-1 0 0]; % A(g)
nu(2,:) = [0 0 1]; % B(g)
nu(3,:) = [1 -1 0]; % A*
nu(4,:) = [0 1 -1]; % B*

% Site balance (unitless)
theta_star = 1 - sum(x(3:4));

% Rate expressions (1/s)
r(1,1) = k(1)*x(1)*Rg*T/Po*theta_star - k(2)*x(3); % A(g) + * <--> A*
r(2,1) = k(3)*x(3) - k(4)*x(4); % A* <--> B*
r(3,1) = k(5)*x(4) - k(6)*x(2)*Rg*T/Po*theta_star; % B* <--> B(g) + *

% Compute xdot:
% Vapor species (M/s)
xdot(1:2,1) = nu(1:2,:)*r*Nsites/V + SV*(Cf' - x(1:2));

% Surface species (1/s)
xdot(3:4,1) = nu(3:4,:)*r;
end

function [SV,tos,TOFbavg,Theta_avg,qsve,PlotData] = square_waveform_CSTR(Nosc,k,T,SV,Cf,tau,x0,C,qdot,UL,UR,BEa0,ssTol)

% previous tsv,xsv,SV,tos,TOFbavg,Theta_avg,BEasv
% ODE solver options
options = odeset('RelTol',1e-6,'AbsTol',1e-8);

% Parameters for data analysis
% Catalyst weight
w = 138; % mg
% Active site density
srho = 20; % micromol/g
% Number of active sites (gmol)
Nsites = w*0.001*srho*0.000001;

% Preallocate matrices
tos = Nosc*mean(tau);
tsv = []; % s
xsv = [];
BEasv = [];
BEadd = [];

z = 1;


% Iterate until convergence
for n = 1:9e18 % can use a while loop instead

    % Preallocate matrices -- for assessing conditions for convergence
    tsvm = 0; % s time vector for middle of run (nos/2 nos/2 - 1)
    xsvm = [];
    tsve = 0; % s time vector or end of run (nosc, nosc-1)
    xsve = [];
    
    % Iterate over oscillations
    % take parameters to do the square waveform oscillation alternation bn
    % L & R enpoint
    for i = 1:Nosc
    
    % Strong binding state
        if bitget(i,1) == 1
    
            % Obtain CSTR data
            [t,x] = ode15s(@(t,x) cstr(t,x,k(1,:),T,SV,Cf),[0 tau(1)],x0(:,i),options);
            
                % Calculate & store BEa-on-stream data
                BEa = BEa0 + UR;
                len = length(t);
                BEadd(1:len,:) = BEa;
                BEasv = [BEasv; BEadd];
                BEadd = []; % Clear after saving
                

                % Store all remaining data
                if isempty(tsv) == 1
                    tsv = t;
                else
                    tsv = [tsv;t + tsv(end)]; % s
                end
                xsv = [xsv;x];


            
            % Store mid & endpoint data for analysis
            if i == round(Nosc/2) || i == round(Nosc/2) - 1
                % Store midpoint data
                tsvm = [tsvm;t + tsvm(end)]; % s
                xsvm = [xsvm;x];
            else
                % Store endpoint data
                if i == round(Nosc) || i == round(Nosc) - 1
                    tsve = [tsve;t + tsve(end)]; % s
                    xsve = [xsve;x];
                end
            end
            
            % Set initial conditions for next oscillation
            x0(:,i + 1) = x(end,:)';
        
            
    
    % Weak binding state
        else
            
            % Obtain CSTR data
            [t,x] = ode15s(@(t,x) cstr(t,x,k(2,:),T,SV,Cf),[0 tau(2)],x0(:,i),options);
                
                % Calculate & store BEa-on-stream data
                BEa = BEa0 + UL;
                len = length(t);
                BEadd(1:len,:) = BEa;
                BEasv = [BEasv; BEadd]; 
                BEadd = []; % Clear after saving
                       
        
                % Store all remaining data
                tsv = [tsv;t + tsv(end)]; % s
                xsv = [xsv;x];
            
            % Set initial conditions for next oscillation
            x0(:,i + 1) = x(end,:)';
            
            
            % Store midpoint & endpoint data for analysis
            if i == round(Nosc/2) || i == round(Nosc/2) - 1
                % Store midpoint data
                tsvm = [tsvm;t + tsvm(end)]; % s
                xsvm = [xsvm;x];
            else
                % Store endpoint data
                if i == round(Nosc) || i == round(Nosc) - 1
                    tsve = [tsve;t + tsve(end)]; % s
                    xsve = [xsve;x];
                end
            end
           
        end
    end

    
   
% Data analysis:
    % Compute conversion (%):
    Xam = (Cf(1) - xsvm(:,1))/Cf(1)*100; % Middle
    Xae = (Cf(1) - xsve(:,1))/Cf(1)*100; % End
    
    % Remove extraneous time point
    tsvm = tsvm(2:end); % s
    tsve = tsve(2:end); % s
    
    % Compute time-averaged conversion (%):
    Xa_avgm = trapz(tsvm,Xam)/range(tsvm);
    Xa_avge = trapz(tsve,Xae)/range(tsve);
    
    % Check steady state
    if abs(Xa_avgm - Xa_avge) > ssTol % 0.0001 loosest tolerance recommended
        % Modify initial conditions
        x0(:,1) = x(end,:)';
        % Track added time (s)
        tos = tos + t(end)*Nosc;
    else
        % Check target conversion
        if abs(Xa_avge - C) > 0.01
            % Modify parameters
            SV = SV*abs(Xa_avge)/C; % 1/s
            qdot = qdot*abs(Xa_avge)/C; % L/s
            % Reset matrices
            x0(:,1) = [Cf';0;0];
            tos = Nosc*mean(tau);
            tsv = []; % s
            xsv = [];
            BEasv = [];
        else
            % Exit infinite loop after 5 iterations at steady-state
            qsve = qdot;
            z = z + 1; % Count how long at steady-state
            if z == 5
                break
            else
            end

        end
    end
end


% Data analysis:
    % Compute production rate (1/s):
    TOFbe = (xsve(:,2) - Cf(2))*qdot/Nsites;

    % Compute time-averaged production rate (1/s):
    TOFbavg = trapz(tsve,TOFbe)/range(tsve);
    % *** these use xsve b/c only want values once reached steady state!

    % Compute time-averaged coverage (unitless):
    Theta_avg(1) = trapz(tsve,xsve(:,3))/range(tsve); % A*
    Theta_avg(2) = trapz(tsve,xsve(:,4))/range(tsve); % B*
    Theta_avg(3) = trapz(tsve,1 - xsve(:,3) - xsve(:,4))/range(tsve); % *

    
% Save last few oscillations for creating plots
    len = length(TOFbe)*5;
    TOFb_inst = (xsv(end-len:end,2) - Cf(2))*qsve/Nsites;

    % Put all plot data into matrix
    PlotData = [tsv(end-len:end,:), xsv(end-len:end,:), BEasv(end-len:end,:),TOFb_inst];
end