% Modified version of the MATLAB template file SolvIVDifI.m used in the
% solution of Example 26.1 of "A First Course on Kinetics and Reaction
% Engineering."
%
function Example_26_1
    % Known quantities and constants
    sp = 0.55;
    Dp = 0.25 * 0.0254; % m
    rho = 32 / 0.02832; % lbm/m^3
    eps = 0.65;
    L = 6; % m
    D = 7 / 100; % m
    U = 100 * 251.996/1.8; % cal/(h m^2 K)
    VFR0 = 200 * 0.02832; % m^3/h
    m0 = 65; % lbm/h
    yA0 = 0.15;
    yB0 = 0.15;
    yI0 = 0.7;
    P0 = 5; % atm
    T0 = 450 + 273.15; % K
    Te = 450 + 273.15; % K
    k0 = 2.9E8; % mol/(h lbm atm^0.5)
    E = 30200; % cal/mol
    mu = 0.3 * 2.41909 / 0.3048; % lbm/(m h)
    dHfA = -70000; % cal/mol
    dHfB = 0; % cal/mol
    dHfZ = -90000; % cal/mol
    CpA = 12.9; % cal/(mol K)
    CpB = 6.7; % cal/(mol K)
    CpZ = 15.7; % cal/(mol K)
    CpI = 5.1; % cal/(mol K)
    % Universal Constants
    Re = 1.987; % cal/(mol K)
    Rpv = 0.08206E-3; % m3 atm/(mol K)
    gc = 32.174*2116.2*60^2*60^2/.3048; % lbm/(m h^2 atm)
    % Additional Calculated Constants
    ntot0 = P0*VFR0/Rpv/T0; % mol/h, ideal gas law
    nA0 = yA0*ntot0; % mol/h
    nB0 = yB0*ntot0; % mol/h
    nI0 = yI0*ntot0; % mol/h
    nI = nI0; % mol/h
    nZ0 = 0; % mol/h
    dH298 = 2*dHfZ - 2*dHfA - dHfB; % cal/mol, equation (17)
    G = 4*m0/pi()/D^2; % lbm/(m^2 h), equation (19)

	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        nA = z(1);
        nB = z(2);
        nZ = z(3);
        T = z(4);
        P = z(5);
        PA = nA*P/(nA + nB + nZ + nI); % equation (15)
        PB = nB*P/(nA + nB + nZ + nI); % equation (16)
        r1 = rho*k0*exp(-E/Re/T)*PA^0.5*PB; % equation (2)
        dH = dH298 + (2*CpZ - 2*CpA - CpB)*(T - 298); % equation (18)
        rho = m0*P/Rpv/T/(nA + nB + nZ + nI); % lbm/m^3, equation (20)
        dzdt = [
            -pi()*D^2/2*r1;
            -pi()*D^2/4*r1;
            pi()*D^2/2*r1;
            (pi()*D*U*(Te-T) - pi()*D^2/4*r1*dH)/(nA*CpA + nB*CpB + nZ*CpZ + nI*CpI);
            -(1-eps)/eps^3*G^2/rho/sp/Dp/gc*(150*(1-eps)*mu/sp/Dp/G + 1.75);
        ];
    end % of internal function odeqns
        
    % Initial and final values
    t0 = 0;
    z0 = [
        nA0
        nB0
        nZ0
        T0
        P0
    ];
    tf = L;
    
    % 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,:);
	
    % Calculate and display the requested results
    T_outlet_Kelvin = z(4)
    pct_conv = 100*(nA0-z(1))/nA0
    P_drop_atm = P0 - z(5)
    
end % of Example_26_1.m
