%% -- Simple example on calculation of                -- %%     
%% -- equilibrum composition for one reaction         -- %%
%% -- CH4 + CO2 <=> 2CO + 2H2    Hr(298K) = 247kJ/mol -- %%
clc;clear;
nc   = 4;          % Number of compound
stoi = [-1 -1 2 2];% Soi. matrix
K    = 4.0e-6;     % Eq. constant at Tref
HR   = 247e3;      % Reaction ethalpy [J/mol]
RG   = 8.314;      % Universal gas constant [J/mol/K]
Tref = 1200;       % Reference temperature [K]
p0   = 1.0;        % Pressure [atm]
xi0  = 0.1;        % Initial extimate of reac. extent
n0   = [1 1 0 0];  % Initial composition [mol]
nstep= 10;         % Number of steps (T)
Tstep= 20;         % Temperature incement [K]

% ... Define n(i) as function of xi
n_i = @(x,i)(n0(i) + stoi(i)*x);

% ... Define K as function of temperature (HR independent on T!)
KT = @(T)(K.*exp(-HR./RG.*(1./T-1./Tref)));

% ... Define K as function of xi
nprod = @(x)(prod(n_i(x,1:nc).^(stoi)));
nsum  = @(x)(sum(n_i(x,1:nc)));
funk  = @(x,T,p)(KT(T) - (p./nsum(x)).^(sum(stoi)).*nprod(x));

% ... Solve
fprintf('%9s %9s %6s %9s\n','Temp','Xi_out','Flag','X(1)');
T = Tref;
p = p0;
for i = 1:nstep
    % Since fzero accepts olny one argument of function funk
    % we have to use following workaround to pass two more
    % arguments "T" and "p" inside the funk (@(x) funk(x,T,p) ... 
    [xi_out,fval,flag] = fzero(@(x) funk(x,T,p),xi0);
    if (flag == 1)
        msg = 'OK';
    else
        msg ='FAILED';
    end
    X1 = (n0(1)-n_i(xi_out,1))/n0(1);
    fprintf('%9.0f %9.3e %6s %9.4f\n',T,xi_out,msg,X1);
    xplot(i) = X1;
    Tplot(i) = T;
    T = T + Tstep;
end
plot(Tplot,xplot);

