% MATLAB function used in the solution of Example 32.2 of "A First Course 
% on Kinetics and Reaction Engineering. It calls separate functions that
% evaluate the equilibrium batch reactor conversion, the first stage of the
% semi-batch operation and the second stage of the semi-batch operation and
% uses the results to answer the questions posed in the problem.
%
function Example_32_2
    % Known quantities and constants (in min, mol, L, atm and K)
    K = 8e-5;
    hZ = 8e-4; % mol/L/atm
    k = 5e-4; % L/mol/min
    V = 100*3.785; % L
    CA0 = 4.; % mol/L
    CB0 = 4.; % mol/L
    T = 80 + 273.15; % K
    Vhead = 8*3.785; % L
    Pmax = 1.5; % atm
    % Constants found in handbooks and other reference sources
    R = 0.08206; % L*atm/mol/K
    % Calculated constants
    nA0 = CA0*V;
    nB0 = CB0*V;
    nY0 = 0.0;
    nZ0 = 0.0;
    % Chosen duration for semi-batch processing
    ttot = 600; % min
    
    % Report the batch conversion (and the final batch pressure)
    [fA_batch,P_batch] = Example_32_2_equil;
    display('Equilibrated Batch Reactor')
    pct_conversion = 100*fA_batch
    final_pressure_atm = P_batch
    
    display('Semi-Batch Reactor')
    % Analyze the first stage of semi-batch processing
    [ts1,nAs1,nBs1,nYs1,nZs1] = Example_32_2_semi_a;
    npoints = length(ts1);
    first_stage_duration_min = ts1(npoints)
    pct_conv = (nA0 - nAs1)*100/nA0;
    first_stage_pct_conversion = pct_conv(npoints)
    % Save the final values from stage 1 to use as initial values for stage
    % 2
    ts1f = ts1(npoints);
    nAs1f = nAs1(npoints);
    nBs1f = nBs1(npoints);
    nYs1f = nYs1(npoints);
    nZs1f = nZs1(npoints);
    
    % Analze the second state of semi-batch processing
    [ts2,nAs2,nBs2,nYs2,nZs2] = Example_32_2_semi_b(ts1f,nAs1f,nBs1f,nYs1f,nZs1f,ttot);
    % Combine the results for the two stages
    ts = [
        ts1
        ts2
    ];
    nAs = [
        nAs1
        nAs2
    ];
    pct_conv = (nA0 - nAs)*100/nA0;
    % Make the requested plot
    plot(ts,pct_conv)
    xlabel('Time (min)')
    ylabel('Percent Conversion')

end % of Example_32_2.m

