% Modified version of the MATLAB template file SolvIVDifI.m used in the
% solution of Example 35.1 of "A First Course on Kinetics and Reaction
% Engineering."
%
function z = Example_35_1_PFR(reactor,q_R2)
    % Known quantities and constants
    L_actual = 6; % m
    D_actual = 7 / 100; % m
    T = 450 + 273.15; % K
    P = 5; % atm
    VFR_FS1 = 200 * 0.02832; % m^3/h
    yA_FS1 = 0.15;
    yB_FS1 = 0.15;
    yI_FS1 = 0.7;
    yZ_FS1 = 0.0;
    k1_R1 = 72; % mol/(h m^3 atm^0.5)
    k1_R2 = 59.5; % mol/(h m^3 atm^0.5)
    
    % Universal Constants
    Rpv = 0.08206E-3; % m3 atm/(mol K)
    
    % Additional Calculated Constants
    ntot_FS1 = P*VFR_FS1/Rpv/T; % mol/h, ideal gas law
    nA_FS1 = yA_FS1*ntot_FS1; % mol/h
    nB_FS1 = yB_FS1*ntot_FS1; % mol/h
    nZ_FS1 = yZ_FS1*ntot_FS1; % mol/h
    nI_FS1 = yI_FS1*ntot_FS1; % mol/h
    V_actual = pi()*D_actual^2*L_actual/4; % m3
    V_R2 = 0.05*V_actual; % m3
    V_R1 = 0.95*V_actual; % m3
    
    % Set constants, initial conditions and final condition for the
    % reactor being modeled
    if (reactor == 'R1')
        k1 = k1_R1;
        tf = V_R1;
        ntot_FS2 = (1-q_R2)*ntot_FS1;
        nI = yI_FS1*ntot_FS2;
        z0 = [
            yA_FS1*ntot_FS2;
            yB_FS1*ntot_FS2;
            yZ_FS1*ntot_FS2;
        ];
    elseif (reactor == 'R2')
        k1 = k1_R2;
        tf = V_R2;
        ntot_FS3 = (q_R2)*ntot_FS1;
        nI = yI_FS1*ntot_FS3;
        z0 = [
            yA_FS1*ntot_FS3;
            yB_FS1*ntot_FS3;
            yZ_FS1*ntot_FS3;
        ];
    else
        display(['Unknown Reactor: ' reactor])
    end

	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        nA = z(1);
        nB = z(2);
        nZ = z(3);
        CA = nA*P/((nA + nB + nZ + nI)*Rpv*T);
        CB = nB*P/((nA + nB + nZ + nI)*Rpv*T);
        r = k1*sqrt(CA)*CB;
        if (CA <= 0.0)
            r = 0.0;
        end
        dzdt = [
            -2*r;
            -r;
            2*r;
        ];
    end % of internal function odeqns
        
    % Initial values of independent variable (other initial and final
    % values were specified above
    t0 = 0;
    
    % Solve the ODEs
    [t,zz] = ode45(@odeqns,[t0 tf],z0);

    % Return the values of the independent and dependent variables when the
    % integration stopping condition was satisfied
	last_value = length(t);
	t_f = t(last_value);
	z = zz(last_value,:);
    
end % of Example_35_1_PFR.m
