% Modified version of the AFCoKaRE MATLAB template file SolvIVDifD.m used
% in the solution of Example 32.2 of "A First Course on Kinetics and
% Reaction Engineering. More specifically, this file performs the
% calculations for the first stage of semibatch operation where the
% pressure relief valve is closed.
%
function [ts1,nAs1,nBs1,nYs1,nZs1] = Example_32_2_semi_a
    % Known quantities and constants (in min, mol, L, atm and K)
    K = 8e-5;
    hZ = 8e-4; % mol/L/atm
    k = 5e-4; % L/mol/min
    V = 100*3.785; % L
    CA0 = 4.; % mol/L
    CB0 = 4.; % mol/L
    T = 80 + 273.15; % K
    Vhead = 8*3.785; % L
    Pmax = 1.5; % atm
    % Constants found in handbooks and other reference sources
    R = 0.08206; % L*atm/mol/K
    % Calculated constants
    nA0 = CA0*V;
    nB0 = CB0*V;
    nY0 = 0.0;
    nZ0 = 0.0;

	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        nA = z(1);
        nB = z(2);
        nY = z(3);
        nZ = z(4);
        CA = nA/V;
        CB = nB/V;
        CY = nY/V;
        CZ = nZ/V;
        r = k*CA*CB*(1 - CY*CZ/K/CA/CB);
        nZdot = -V*r/(1 + V*hZ*R*T/Vhead);
        dzdt = [
            -V*r
			-V*r
            V*r
            nZdot + V*r
        ];
    end % of internal function odeqns
        
    % Initial values
    t0 = 0.0;
    z0 = [
        nA0
        nB0
        nY0
        nZ0
    ];

    tf = 100.; % should never reach this
    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;
        
        nZ = z(4);
        CZ = nZ/V;
        PZ = CZ/hZ;
        stop_when = PZ - Pmax;
    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
	
    % Set the return values
    ts1 = t;
    nAs1 = zz(:,1);
    nBs1 = zz(:,2);
    nYs1 = zz(:,3);
    nZs1 = zz(:,4);
    
end % of Example_32_2_semi_a.m
