% Modified version of the AFCoKaRE MATLAB template file SolvIVDifD.m used
% in the solution of Activity 32.1 of "A First Course on Kinetics and
% Reaction Engineering. More specifically, this file performs the
% calculations for the batch reactor.
%
function Activity_32_1_batch
    % Known quantities and constants (in min, mol, L)
    k1 = 8e-3; % L/mol/min
    k2 = 4e-3; % L/mol/min
    Vtot = 2000.0; % L
    CA0 = 3.0; % mol/L
    CB0 = 1.0; % mol/L
    VA = 500.0; % L
    VB = 1500.0; % L
    fB = 0.95;
    % Calculated constants
    nA0 = CA0*VA;
    nB0 = CB0*VB;
    nD0 = 0.0;
    nU0 = 0.0;

	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        nA = z(1);
        nB = z(2);
        nD = z(3);
        nU = z(4);
        CA = nA/Vtot;
        CB = nB/Vtot;
        r1 = k1*CA*CB;
        r2 = k2*CB^2;
        dzdt = [
            -Vtot*r1
            Vtot*(-r1 - 2*r2)
            Vtot*r1
            Vtot*r2
        ];
    end % of internal function odeqns
        
    % Initial values
    t0 = 0.0;
    z0 = [
        nA0
        nB0
        nD0
        nU0
    ];

    % Final value
    tf = 1000.;
    options = odeset('Events',@stop);
	[t, zz, te, ze, ie] = ode45(@odeqns,[t0, tf],z0,options);
    % NOTE: For stiff equations replace ode45 with ode15s
    
    % Function that provides the integration stopping criterion
    function [stop_when, isterminal, direction] = stop(t,z)
        isterminal = 1;
        direction = 0;
        
        stop_when = z(2) - nB0*(1-fB);
    end % of internal function stop

    % Extract the final values of the independent and dependent variables
	last_value = length(t);
	t_f = t(last_value); % final value of the independent variable
	z = zz(last_value,:); % final values of the dependent variables
    
    % Warn if the independent variable reached tf above
    if (t_f >= tf)
        display('WARNING: the final value of the independent variable was')
        display('reached, instead of the final value of the dependent')
        display('variable.');
        display('Set tf to a larger value and re-execute the function')
    end
	
    % Calculate the selectivity
    sel_D_per_U_batch = z(3)/z(4)
    t_batch = t_f
    
end % of Activity_32_1_batch.m
