% 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 semi-batch reactor.
%
function Activity_32_1_semi
    % Known quantities and constants (in min, mol, L)
    k1 = 8e-3; % L/mol/min
    k2 = 4e-3; % L/mol/min
    CA0 = 3.0; % mol/L
    CB0 = 1.0; % mol/L
    V0 = 500.0; % L
    Vmax = 2000.0; % L
    VFR = 150.0/60; % L/min
    fB = 0.95;
    % Calculated constants
    nA0 = CA0*V0;
    nB0 = 0.0;
    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);
        V = V0 + VFR*t;
        nBdot = VFR*CB0;
        if V > Vmax
            V = Vmax;
            nBdot = 0.0;
        end
        CA = nA/V;
        CB = nB/V;
        r1 = k1*CA*CB;
        r2 = k2*CB^2;
        dzdt = [
            -V*r1
            nBdot + V*(-r1 - 2*r2)
            V*r1
            V*r2
        ];
    end % of internal function odeqns
        
    % Initial values
    t0 = 0.0;
    z0 = [
        nA0
        nB0
        nD0
        nU0
    ];

    % Final value
    tf = 2000.;
    options = odeset('Events',@stop);
	[t, zz, te, ze, ie] = ode15s(@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;
        
        nBadded = CB0*(VFR*t);
        V = V0 + VFR*t;
        if V > Vmax
            nBadded = CB0*(Vmax-V0);
        end
        stop_when = z(2) - nBadded*(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. Set tf to a larger value and re-execute')
        display('    the function')
    end
	
    % Calculate the selectivity
    sel_D_per_U_batch = z(3)/z(4)
    t_semi = t_f
    
end % of Activity_32_1_semi.m
