% Modified version of the AFCoKaRE MATLAB template file SolvNonDif.m used
% in the solution of Example 32.2 of "A First Course on Kinetics and
% Reaction Engineering. More specifically, this file performs the
% equilibrium calculations for the semi-batch reactor without a relief 
% valve.
%
function [fA_equil,Pfinal] = Example_32_2_equil
    % 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 equations
    function f = evalEqns(z)
        CA = z(1);
        CZ = z(2);
        nA = CA*V;
        extent = nA0 - nA;
        nB = nB0 - extent;
        nY = nY0 + extent;
        CB = nB/V;
        CY = nY/V;
        nZtot = nZ0 + extent;
        nZliq = CZ*V;
        nZhead = nZtot - nZliq;
        PZ = nZhead*R*T/Vhead;
        f = [
            K*CA*CB - CY*CZ
            PZ*hZ - CZ
        ];
    end % of internal function evalEqns
    
    % guesses for the solution
    z_guess = [
        0.9*CA0
        0.1*CA0
    ];

    % Solve the set of algebraic equations
    z = fsolve(@evalEqns, z_guess);
    display('The solver found the following values for the unknowns:');
    z
    display('The corresponding values of the functions being solved are as follows:');
    f = evalEqns(z)
    
    % calculate the conversion and set the return variable values
    CA = z(1);
    CZ = z(2);
    nA = CA*V;
    extent = nA0 - nA;
    nZtot = nZ0 + extent;
    nZliq = CZ*V;
    nZhead = nZtot - nZliq;
    fA_equil = (nA0 - nA)/nA0;
    Pfinal = nZhead*R*T/Vhead;
	
end % of Example_32_2_equil.m