% CSTR_2D_NAE.m
% 2-D example of CSTR calculation with homotopy
% K. Beers. MIT ChE. 10.34. 7/14/03

function iflag_main = CSTR_2D_NAE();

iflag_main = 0;

% First, set values of system parameters
Param.Q = 0.1;  Param.V = 1;
Param.A_in = 1;  Param.B_in = 2;
Param.C_in = 0;  Param.D_in = 0;
Param.k1 = 1;  Param.k2 = 1;
% Ask user whether to use homotopy
iuse_homotopy = input('Use homotopy 0=no, 1=yes : ');
% set options for fsolve() MATLAB solver
Options = optimset('LargeScale', 'off');
% make initial guess using inlet conc.
x0 = [Param.A_in; Param.B_in; Param.C_in; Param.D_in];

if(~iuse_homotopy)  % solve without homotopy
    [x,fval,exitflag] = fsolve(@CSTR_2D_NAE_calc_f,x0,Options,Param);
        
else  % use homotopy
    Q_set = Param.Q;  Q_max = Param.Q * 1000;  num_homotopy = 5;
    Q_vect = logspace(log10(Q_max),log10(Q_set),num_homotopy);
    for iter=1:num_homotopy
        Param.Q = Q_vect(iter);
        [x,fval,exitflag] = fsolve(@CSTR_2D_NAE_calc_f,x0,Options,Param);
        if(exitflag<=0)
            disp(['Not coverged at iter = ', int2str(iter)]);
            return;
        end
        x0 = x;  % use last result as next initial guess
    end
end

% display result
if(exitflag==0)
    disp('Max. # of iterations reached without convergence');
    return;
elseif(exitflag<0)
    disp(['Error with exitflag = ', int2str(exitflag)]);
    return;
else
    disp('Converged to a solution');  disp(' ');
    disp('The steady state concentrations are : ');
    disp(['[A] = ', num2str(x(1))]);
    disp(['[B] = ', num2str(x(2))]);
    disp(['[C] = ', num2str(x(3))]);
    disp(['[D] = ', num2str(x(4))]);
end

iflag_main = 1;
return;


%=======================
function f = CSTR_2D_NAE_calc_f(x,Param);

f = zeros(4,1);
f(1) = Param.Q*(Param.A_in - x(1)) + ...
    Param.V*(-Param.k1*x(1)*x(2));
f(2) = Param.Q*(Param.B_in - x(2)) + ...
    Param.V*(-Param.k1*x(1)*x(2) - Param.k2*x(3)*x(2));
f(3) = Param.Q*(Param.C_in - x(3)) + ...
    Param.V*(Param.k1*x(1)*x(2) - Param.k2*x(3)*x(2));
f(4) = Param.Q*(Param.D_in - x(4)) + ...
    Param.V*(Param.k2*x(3)*x(2));

return;