% Modified version of SolvIVDifI used to model the PFR in Activity 30.1 
% from "A First Course on Kinetics and Reaction Engineering."
%
function [nAout,nWout,nPout,Tout] = Activity_30_1_pfr(Tin)
    % 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)
    nAin = 0.07; % kmol/s
    nWin = 1.67; % kmol/s
    nPin = 0; % kmol/s
    VFR = 0.04; % m^3/s

	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        nA = z(1);
        nW = z(2);
        nP = z(3);
        T = z(4);
        CA = nA/VFR;
        CW = nW/VFR;
        CP = nP/VFR;
        k = k0*exp(-EoverR/T);
        K = K0*exp(HoverR/T);
        r = k*CA*CW*(1-CP/K/CA/CW);
        dH = dH298 + (CpP - CpA - CpW)*(T - 298);
        dzdt = [
            -r
            -r
            r
            -r*dH/(nA*CpA + nW*CpW + nP*CpP)
        ];
    end % of internal function odeqns
        
    % Initial and final values
    t0 = 0;
    z0 = [
        nAin
        nWin
        nPin
        Tin
    ];
    tf = V;

    % Solve the ODEs
    [t,zz] = ode45(@odeqns,[t0 tf],z0);

    % Extract the final values of the independent and dependent variables
    last_value = length(t);
    t_f = t(last_value);
    z = zz(last_value,:);

    nAout = z(1);
    nWout = z(2);
    nPout = z(3);
    Tout = z(4);
        
end % of Activity_30_1_pfr.m
