% Modified version of the MATLAB template file SolvIVDifD.m used in the
% solution of Example 20.1 of "A First Course on Kinetics and Reaction
% Engineering." This function models the first step of the operational
% protocol described in that example.
%
function [t_f,z] = Example_20_1_a
    % Known quantities and constants (cal, mol, K, min, ft2, kg, L)
    k0 = 2.59e9; % min-1
    E = 16500; % cal/mol
    R = 1.987; % cal/mol/K
    dH = -22200; % cal/mol
    V = 4; % L
    CA0 = 2; % mol/l
    nA0 = CA0*V; % mol
    Te_coil = 120 + 273.15; % K
    Te0_jacket = 23 + 273.15; % K
    cp = 440; % cal/L/K
    Ve = 0.5; %L
    A_jacket = 0.6; % ft2
    U_jacket = 1.13e4/60; % cal/ft2/min/K
    A_coil = 0.23; % ft2
    U_coil = 3.8e4/60; % cal/ft2/min/K
    rho = 1; % kg/L
    cpe = 1000; % cal/kg/K
    T0 = 23 + 273.15; % K
    Tf = 50 + 273.15; % K
    m = 0; % kg/min (no cooling water flowing)

	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        nA = z(1);
        nB = z(2);
        T = z(3);
        Te_jacket = z(4);
        r1 = k0*exp(-E/R/T)*nA/V;
        Q = U_jacket*A_jacket*(Te_jacket-T)+U_coil*A_coil*(Te_coil-T);
        dzdt = [
            -r1*V
            r1*V
            (Q-V*dH*r1)/V/cp
            (m*cpe*(Te0_jacket-Te_jacket)-U_jacket*A_jacket*(Te0_jacket-T))/rho/Ve/cpe
        ];
    end % of internal function odeqns
        
    % Initial values
    t0 = 0.0;
    z0 = [
        nA0
        0
        T0
        Te0_jacket
    ];
    tf = 5000.0;
    options = odeset('Events',@stop);
	[t, zz, te, ze, ie] = ode45(@odeqns,[t0, tf],z0,options);
    
    % Function that provides the integration stopping criterion
    function [stop_when, isterminal, direction] = stop(t,z)
        isterminal = 1;
        direction = 0;
        
        stop_when = Tf - z(3);
    end % of internal function stop

    % 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,:);
	
    % Plot CA vs. t
    CA = zz(:,1)/V;
    figure
    plot(t,CA)
    title(['Heating Step'])
    xlabel('t (min)')
    ylabel('CA (M)')
    % Plot Te and T vs. t
    figure
    plot(t,zz(:,3),t,zz(:,4))
    title(['Heating Step'])
    xlabel('t (min)')
    ylabel('T (K)')
    legend('Reactor','Jacket')
    
end % of Example_20_1_a
