% Modified version of the AFCoKaRE MATLAB template file SolvIVDifD.m used
% in the solution of Example 32.1 of "A First Course on Kinetics and
% Reaction Engineering.
%
function Example_32_1
    % Known quantities and constants (in min, mol, J, K, cc, g)
    k0 = 1.192E15; % /min
    E = 97600.; % J/mol
    dH = -58615.; % J/mol
    Te = 60 + 273.15; % K
    UA = 260*4.184; % J/min/K
    VB0 = 67.0; % cc
    Vcat = 0.3; % cc
    VZ0 = 283.0; % cc
    T0 = 60 + 273.15; % K
    Cp_soln = 2.68; % J/cc/K
    Tstream = 21 + 273.15; % K
    VFR = 64.5; % cc/min
    t_semi = 5.5; % min
    Tfinal = 60 + 273.15; % K
    P = 1.0; % atm
    % Data from handbooks, etc.
    R = 8.31446; % J/mol/K
    rhoA = 1.082; % g/cc
    MWA = 102.; % g/mol
    CpA = 168.2; % J/mol/K
    rhoB = 1.0; % g/cc
    MWB = 18.; % g/mol
    rhoZ = 1.049; % g/cc
    MWZ = 60.; % g/mol
    % Calculated constants
    V0 = VB0 + VZ0 + Vcat;
    nA0 = 0;
    nB0 = VB0*rhoB/MWB;
    nZ0 = VZ0*rhoZ/MWZ;

	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        % t is the value of the independent variable
        % z is a vector containing the values of the dependent variables
        nA = z(1);
        nB = z(2);
        nZ = z(3);
        T = z(4);
        if (t > t_semi)
            ndotA = 0;
            V = V0 + VFR*t_semi;
            dVdt = 0;
        else
            ndotA = VFR*rhoA/MWA;
            V = V0 + VFR*t;
            dVdt = VFR;
        end
        CA = nA/V;
        r = k0*exp(-E/R/T)*CA;
        rxnHeat = -V*r*dH;
        sensHeat = ndotA*CpA*(T-Tstream);
        Qdot = UA*(Te-T);
        expansWork = P*dVdt*0.1013; % J/min
        netHeat = rxnHeat + expansWork - sensHeat + Qdot;
        dzdt = [
            ndotA - V*r;
            -V*r;
            2*V*r;
            (-V*r*dH + expansWork - sensHeat + Qdot)/(V*Cp_soln);
        ];
    end % of internal function odeqns
        
    % Initial values
    t0 = 0;
    z0 = [
        nA0
        nB0
        nZ0
        T0
    ];
    tf = 10.0;
    
    % Solve the ODEs
    [t, zz] = ode15s(@odeqns,[t0, tf],z0);
    
    % Extract the final values of the independent and dependent variables
	last_value = length(t);
	t_f = t(last_value); % final value of the independent variable
	z = zz(last_value,:); % final values of the dependent variables
    
    % calculate the concentration versus time
    conc_A = zeros(1,length(t));
	for i=1:length(t);
        if t(i) < t_semi
            V = V0 + VFR*t(i);
        else
            V = V0 + VFR*t_semi;
        end
        conc_A(i) = zz(i,1)/V;
    end
    % convert to degrees C and mol/L
    T_C = zz(:,4)-273.15;
    conc_A = conc_A*1000.0;
    % make the plots
    figure
    plot(t,T_C)
    xlabel('Time (min)')
    ylabel('T (deg C)')
    figure
    plot(t,conc_A)
    xlabel('Time (min)')
    ylabel('Concentration of A (mol/L)')
    
    % calculate the batch-wise initial concentration of A (mol/L)
    nA_batch = VFR*t_semi*rhoA/MWA;
    V0_batch = V0 + VFR*t_semi;
    CA_batch = nA_batch/V0_batch*1000
    
end % of Example_32_1.m
