% Modified version of the MATLAB template file SolvNonDif.m used for the
% solution of Activity 30.1 from "A First Course on Kinetics and Reaction
% Engineering."
%
function Activity_30_1
    % Known quantities and constants (in consistent units)
    V = 0.5; % m^3
    k0 = 1.2E12; % m3/(kmol s)
    EoverR = 13000; % K
    K0 = 4.2E-15; % m3/kmol
    HoverR = 11300; % K
    dH298 = -86000; % kJ/kmol
    CpA = 412; % kJ/(kmol K)
    CpW = 76; % kJ/(kmol K)
    CpP = 512; % kJ/(kmol K)
    nAa = 0.07; % kmol/s
    nWa = 1.67; % kmol/s
    nPa = 0; % kmol/s
    VFR = 0.04; % m^3/s
    Ta = 300; % K
    
    % Function that evaluates the equations
    function f = evalEqns(z)
        Tb = z(1);
        % Solve the pfr design equations to get stream c properties
        [nAc,nWc,nPc,Tc] = Activity_30_1_pfr(Tb);
        Td = Ta + 5;
        term1 = (nAc*CpA + nWc*CpW + nPc*CpP)*(Tc-Td);
        term2 = (nAa*CpA + nWa*CpW + nPa*CpP)*(Tb-Ta);
        f = [
           term1-term2
        ];
    end % of internal function evalEqns
    
    % guesses for the solution
    z_guess = [
        Ta + 50
    ];

    % 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)
    
    Tb = z(1);
    [nAc,nWc,nPc,Tc] = Activity_30_1_pfr(Tb);
    display('With thermal backmix:')
    outlet_temperature = Ta + 5
    pct_conv = 100*(nAa - nAc)/nAa
    
	% Case with no backmix, feed to reactor is at Ta
    [nAb,nWb,nPb,Tb] = Activity_30_1_pfr(Ta);
    display('Without thermal backmix:')
    outlet_temperature = Tb
    percent_conversion = 100*(nAa-nAb)/nAa

end % of Activity_30_1.m