% MATLAB file used to solve Example 3.1 from "A First Course on Kinetics
% and Reaction Engineering."
function Example_3_1
    % Data given in the problem statement, converted to units of J, mol and
    % K
    r_gas = 8.3144;
    % stoichiometric coefficients (1 = CO, 2 = H2 and 3 = MeOH)
    nu = [-1
        -2
        1];
    % heats of formation at 298 K
    hf = [-110530
    	0
    	-201000];
    % free energies of formation at 298 K
    gf = [-137368.9
        0
        -162615.3];
    % heat capacity constants
    A_co = 25.567590;
    A_h2 = 33.1078;
    A_meoh = 21.138073;
    B_co = 6.096130;
    B_h2 = -11.508;
    B_meoh = 70.878654;
    C_co = 4.054656;
    C_h2 = 11.6093;
    C_meoh = 25.853554;
    D_co = -2.671301;
    D_h2 = -2.8444;
    D_meoh = -28.497905;
    E_co = 0.131021;
    E_h2 = -0.159665;
    E_meoh = 0;
    
    % Expressions for the heat capacities as a cell array of 
    % anonymous functions; x is the temperature in K
    Cp = {@(x) (A_co + B_co*x/1000 + C_co*(x/1000).^2 + D_co*(x/1000).^3 + E_co./(x/1000).^2);
        @(x) (A_h2 + B_h2*(x/1000) + C_h2*(x/1000).^2 + D_h2*(x/1000).^3 + E_h2./(x/1000).^2);
        @(x) (A_meoh + B_meoh*(x/1000) + C_meoh*(x/1000).^2 + D_meoh*(x/1000).^3 + E_meoh./(x/1000).^2);
        };
    
    % Calculate the standard heat and free energy of reaction at 298 K
    h_298 = transpose(nu)*hf;
    g_298 = transpose(nu)*gf;
    
    % Calculate the equilibrium constant at 298 K
    K_298 = exp(-g_298/r_gas/298.);
    
    % Calculate the equilibrium constant at the three specified
    % temperatures using the two internal functions defined below to
    % calculate the heat of reaction at temperature T and the equilibrium
    % constant at temperature T
    K_500 = K_eq(500)
    K_550 = K_eq(550)
    K_600 = K_eq(600)
    
    % Internal function that numerically calculates the heats of reaction
    % at temperatures in the vector T
    function h_T = hr(T)
        h_T = zeros(1,length(T));
        for j=1:length(T)
            % Sensible heat
            h_s = 0.;
            for i=1:length(nu)
                h_s = h_s + nu(i)*integral(Cp{i},298.,T(j));
            end

            % Standard heat of reaction at T in J/mol
            h_T(j) = h_298 + h_s;
        end
    end
    
    % Internal function that numerically calculates the equilibrium
    % constant at temperature T
    function K_T = K_eq(T)
        fun = @(x) (hr(x)/r_gas./x./x);
        integr = integral(fun,298.,T);
        K_T = K_298*exp(integr);
    end
end

