% Modified version of the MATLAB template file SolvIVDifI.m used in the
% solution of Example 26.2 of "A First Course on Kinetics and Reaction
% Engineering."
%
function Example_26_2
    % Known quantities and constants
    D = 10; % cm
    L = 500; % cm
    CA0 = 1.0/1000.0; % mol/cm^3
    CB0 = 1.2/1000.0; % mol/cm^3
    nD0 = 0;
    nS0 = 0;
    nU0 = 0;
    T0 = 30 + 273.15; % K
    k01 = 8.93E2*1000.; % cm^3/(mol min)
    E1 = 6700; % cal/mol
    k02 = 9.24E2*1000.; % cm^3/(mol min)
    E2 = 5100; % cal/mol
    dH1 = -11300; % cal/mol
    dH2 = -11700; % cal/mol
    Cp = 1.0; % cal/(g K)
    rho = 1.0; % g/cm^3
    % Universal Constants
    R = 1.987; % cal/(mol K)
    
    % Create 100 element vectors to hold the data for later plotting and
    % set the range of values of VFR to be used
    VFR = linspace(0,1.0E4);
    concD = zeros(100,1);
    % Create variables to hold the maximum CD and corresponding VFR
    VFR_at_max = 0;
    concD_max = 0;
    % Solve the design equations for the chosen range of VFRs
    for i = 1:100
        % Additional Calculated Constants
        nA0 = VFR(i)*CA0;
        nB0 = VFR(i)*CB0;
        
        % Initial and final values
        t0 = 0;
        z0 = [
            nA0
            nB0
            nD0
            nS0
            nU0
            T0
        ];
        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 the concentration of D
        concD(i) = z(3)./VFR(i);
        if (concD(i) > concD_max)
            concD_max = concD(i);
            VFR_at_max = VFR(i);
        end
    end % of loop through range of values of VFR
    
    % Generate the plot
    figure
    plot(VFR,concD)
    xlabel('Vol. Flow Rate (cm^3/min)')
    ylabel('Concentration of D (mol cm^3)')
    
    % Report the VFR and CD_max
    concD_max
    VFR_at_max
    
    % Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        nA = z(1);
        nB = z(2);
        nD = z(3);
        nS = z(4);
        nU = z(5);
        T = z(6);
        CA = nA/VFR(i);
        CB = nB/VFR(i);
        CD = nD/VFR(i);
        r1 = k01*exp(-E1/R/T)*CA*CB;
        r2 = k02*exp(-E2/R/T)*CA*CD;
        dzdt = [
            pi()*D^2/4*(-r1-r2)
            pi()*D^2/4*(-r1)
            pi()*D^2/4*(r1-r2)
            pi()*D^2/4*(r1+r2)
            pi()*D^2/4*(r2)
            -pi()*D^2/4*(r1*dH1+r2*dH2)/VFR(i)/rho/Cp
        ];
    end % of internal function odeqns

end % of Example_26_2.m
