% Simulates dynamic batch reactor for A-to-B reaction
    % NOTE: need to carefully check that these reach steady-state; tolerance
    % might not be tight enough.
clc
clear


% Specify parameters
    % Scaling parameters
    gamma = -4.5; % Linear scaling slope, (-)
    delta = 1.4; % Equivalent surface enthalpy, (eV)
    alpha = 0.6; % BEP slope, (-)
    beta = 100; % BEP offset, (kJ/mol)

% Waveform parameters
    fosc = 1e-1; % Frequency, (Hz)
    delU = 0.4; % Amplitude, (eV)
    BEa0 = 1.40; % Waveform centerpoint, (eV)

 
% Simulate dynamic batch reactor
	[tos,A_ss,B_ss,Theta_avg,PlotData] = dynamic_catalysis_BR(alpha,beta,gamma,delta,fosc,delU);
    results = [delU,fosc,A_ss,B_ss,Theta_avg]
    
    % Variable names
    %     A_ss, B_ss = time-averaged, dynamic steady-state concentration of A and B (M)
    %     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)
    



function [tos,A_ss,B_ss,Theta_avg,PlotData] = dynamic_catalysis_BR(alpha,beta,gamma,delta,fosc,delU)
% Set constants
% Gas constant
R = 8.31446261815324; % J/mol-K
Rg = R*0.01; % L-bar/mol-K

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

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


% Initial binding energy 
BEa0 = 1.4; % eV

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

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

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

% Initial conditions
t0 = 0;
x0 = [Cf';0;0];


% 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);

% Oscillation endpoints (eV)
UR = delU/2; % Strong binding
UL = UR - delU; % Weak binding


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

% Run dynamic simulation
    tic
    [A_ss,B_ss,Theta_avg,tos,tsv,xsv,BEasv] = square_waveform_batch(Nosc,k,T,Cf,tau,x0,t0,UL,UR,BEa0);
    toc
    
    PlotData = [tsv,xsv,BEasv];
end

function k = batchreactor_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 = batchreactor(t,x,k,T,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;

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

function [A_ss,B_ss,Theta_avg,tos,tsv,xsv,BEasv] = square_waveform_batch(Nosc,k,T,Cf,tau,x0,t0,UL,UR,BEa0)

% ODE solver options
% Strict tolerances recommended for these BR simulations. Loosen these too much
%   and you run the risk of LARGE errors, e.g., rxn directionality wrong.
options = odeset('RelTol',1e-8,'AbsTol',1e-10);


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



% Iterate until convergence
for n = 1:9e18

% 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 = [];
    
% Simulate 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 % i is odd number
    
            % Obtain batchreactor data
            [t,x] = ode15s(@(t,x) batchreactor(t,x,k(1,:),T,Cf),[0 tau(1)],x0(:,i),options);
       
            % Store all data
            if isempty(tsv) == 1
                tsv = t;
            else
                tsv = [tsv;t + tsv(end)]; % s
            end
            xsv = [xsv;x];
            len = length(t);
            BEadd(1:len,:) = BEa0 + UR;            
            BEasv = [BEasv; BEadd];
            BEadd = [];

            % Store midpoint data
            if i == round(Nosc/2) || i == round(Nosc/2) - 1 
                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 batchreactor data
            [t,x] = ode15s(@(t,x) batchreactor(t,x,k(2,:),T,Cf),[0 tau(1)],x0(:,i),options);

            % Store all data
            tsv = [tsv;t + tsv(end)]; % s
            xsv = [xsv;x];
            len1 = length(t);
            BEadd(1:len1,:) = BEa0 + UL;
            BEasv = [BEasv; BEadd];
            BEadd = [];
            
            % Store midpoint data
            if i == round(Nosc/2) || i == round(Nosc/2) - 1
                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,:)';
        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
    mean_Xa = (Xa_avgm + Xa_avge)/2;
    if 100*abs(Xa_avgm - Xa_avge)/mean_Xa > 0.0001  % ss tolerance in relative % difference
        % Modify initial conditions
        x0(:,1) = x(end,:)';
        clear xsve xsvm tsve tsvm
    else
        break
    end
    
end


% Calculate correct time on stream based on IC
    tos = tsv(end) + t0;

% 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;

    % Get steady-state concentrations (M)
    A_ss = trapz(tsve,xsve(:,1))/range(tsve);
    B_ss = trapz(tsve,xsve(:,2))/range(tsve);

    % 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); % *
    
end